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ABSTRACT 


The collection of signals intelligence via passive direction finding and geolocation of radio 
frequency signals is of great concern to the military for its contribution to the development 
of battlespace awareness. Basic subspace direction finding techniques provide a method 
of determining the direction-of-arrival (DOA) of multiple signals on an array of receivers, 
but they have an inherent limitation in that they are narrowband by design. 

The impact of various signal frequencies, bandwidths, and signal to noise ratios 
present in the source signals received by a sparse array using the multiple signals 
classification (MUSIC) subspace direction-finding algorithm are evaluated in this thesis. 
Additionally, two performance enhancements are presented: one that reduces the MUSIC 
computational load and one that provides a method of utilizing collector motion to resolve 
DOA ambiguities. 
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Executive Summary 


Background 

Joint Publication 1-02 defines signals intelligence (SIGINT) as intelligence derived from 
communications, electronic, and foreign instrumentation signals [1]. Direction finding and 
geolocation are important parts of that intelligence; passive direction finding is the process 
of estimating an emitter’s direction from the receiver’s position, and geolocation is then the 
process of using one or more direction-of-arrival (DOA) results in combination with known 
information about receiver or emitter motion and relative geometry in order to determine 
the emitter’s position in terms of an earth-fixed coordinate system. For the military, obtain¬ 
ing location information for emitters is a vital part of understanding an adversary’s order 
of battle and force strength. 

In the future, the signal environment will only become more dense, and the signals 
themselves will become lower power and wider bandwidth. These trends mean that sig¬ 
nals will become more difficult to detect and more likely to overlap in time and bandwidth. 
Airborne and space-based collectors are ideal for conducting SIGINT because they provide 
more area coverage and their motion enables finer location resolution than terrestrial sys¬ 
tems, but they have limitations on the size and number of receivers which can be carried. 
The large coverage area exacerbates the possibility of signal overlap, and the distance can 
make detection of low power signals difficult. Subspace algorithms provide an ability to 
resolve multiple overlapping signals, but many are based on a narrowband assumption or 
require strict receiver antenna array arrangements, which can limit their usefulness [2], [3]. 

Objective 

In order to evaluate its performance with wideband signals and sparse arrays, the multiple 
signals classification (MUSIC) algorithm of [2] has been implemented and examined. 
MUSIC was one of the first subspace algorithms to be developed; many modifications 
and improvements to it have been made since its introduction [4]-[6], but an analysis of 
the performance of the basic algorithm will provide a baseline reference to which other 
methods can be compared. Additionally, it is a goal of this thesis to evaluate how wideband 
the input signals can be before the algorithm breaks down. Two enhancements to the al- 
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gorithm have also been developed: a method of improving the eomputational performanee 
of MUSIC without saerifieing resolution, and a method of utilizing eolleetor motion to 
determine emitter loeation for an array with ambiguities. 


Performance Analysis 

For the majority of the analysis, a sparse array with six elements, a large baseline, and no 
ambiguities was utilized for the eolleetor. For the emitters, variations of three basie signal 
types were modeled: pure sinusoid, linear frequeney modulation (LFM), and quadrature 
phase-shift keying (QPSK). The evaluation was divided into four areas examining the 
effeets on the MUSIC results: signal eenter frequeney, bandwidth, DOA, and signal-to- 
noise ratio (SNR). 

Signals with eenter frequeney other than the eolleetor’s eenter frequeney were found to 
have a shifted DOA in the MUSIC results. This shift ean easily be eorreeted, however, 
assuming that an aeeurate estimation ean be made of the souree signal eenter frequeneies. 
A eorreeted DOA ean be determined using the ratio of the eolleetor eenter frequeney to the 
estimated souree signal’s eenter frequeney. 

Signal bandwidth leads to more eomplieated distortions that eannot be as easily eor¬ 
reeted as ean signal frequeney offsets. The first effeet to note is that the amplitude of the 
MUSIC response deereases with inereased bandwidth, though even at bandwidths up to 40 
pereent of the eenter frequeney a peak of suffieient amplitude above the baekground ean be 
resolved. In testing signals with two widely separated DOAs but signifieantly overlap in 
bandwidth, it was found that MUSIC was easily able to resolve the two signals. A third test 
was eondueted to determine the performanee for wideband signals with elose DOAs but 
separation in frequeney, and MUSIC was again able to resolve the DOAs of two signals. 
In this test, however, the signals were imperfeetly reeonstrueted with blending of the two 
signals present whieh beeame worse with an inerease in bandwidth. 

A final bandwidth analysis with a third signal whieh was separated in DOA but over¬ 
lapped in frequeney with one of the prior two was then evaluated. It was found that when 
the overlap is slight (narrow bandwidth third signal), the DOAs are all able to be resolved 
via MUSIC, but when the overlap is large (wide bandwidth third signal), the MUSIC re¬ 
sponse beeomes degraded: the two elosely spaeed signals ean no longer be resolved and 
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the third signal response can be distorted. 

The effect of DOA on the MUSIC response is related to the bandwidth. For sinusoidal 
signals, the DOA has little to no impact on the MUSIC response, but for wideband signals 
the amplitude of the MUSIC peak can widen and decrease drastically as the offset from 
array boresight increases. 

The last signal parameter evaluated was the effect of SNR. One might expect that the 
MUSIC peak amplitudes would decrease with decreased SNR, and this was found to be the 
case. An examination of the statistics of the Monte Carlo simulation reveals more infor¬ 
mation, at low SNRs the mean estimated DOA is still accurate, but the standard deviation 
increases drastically. A second evaluation of SNR effects was conducted with two closely 
spaced signals, and it was found that the signal DOAs were easily resolvable at 10 dB SNR 
and just resolvable at 0 dB SNR, but at -10 dB SNR there was only a single peak present. 
A final analysis examined the impact of multiple signals with different SNRs. Here band¬ 
width played a role again in leading to additional distortions beyond those of narrowband 
signals. In general, sinusoidal signals could be still be resolved up to a 20 dB SNR differ¬ 
ence, but signals with 10% bandwidth could only be resolved up to a 10 dB difference in 
SNR. 

Computational Enhancement 

One of the major limitations of the MUSIC algorithm is that the search over the angular 
visible range is computationally expensive. The primary driver of that computation is not 
the MUSIC calculation itself, but the fact that the MUSIC calculation must be performed 
over a discrete range of possible signal DOAs. In order to obtain a high resolution estimate 
of the DOA a small angular step size must be used which means a large number of discrete 
angles must be evaluated. A method of reducing the number of discrete angles without 
also reducing the resolution would be to run the MUSIC algorithm more than once. On the 
first run, a large step size should be used which provides only a rough estimate of signal 
direction. Those estimates can then be used to re-run MUSIC with a much smaller step 
size over only the angles surrounding each of the signal direction estimates, with the result 
of obtaining the high resolution DOA estimate desired. With two signals present, it was 
found that a ~ 30x improvement could be obtained on an azimuth-only evaluation with a 



final resolution of approximately 100 ft at a distance of nearly 200 nmi. 


Ambiguity Resolution 

Array ambiguities, which arise when the minimum receive element spacing is too large, 
are unresolvable for fixed arrays using subspace methods. Ambiguity resolution can be 
accomplished, however, when the collector is moving at a known velocity and multiple 
collections of the same signal can be obtained. This is because the true DOA will always 
point in the direction of the emitter, but any ambiguous DOAs will tend to drift. By track¬ 
ing the intersection points of DOA vectors from subsequent collections, it should easily 
be determinable that one set of crossing points remains relatively fixed while another set 
drifts. Based on this, a metric was developed which uses the motion of subsequent DOA 
intersection points and can provide an estimate of the correct DOA in only three collections. 

It was found that this metric was less accurate for wideband and low SNR signals be¬ 
cause in those cases the estimated DOAs of each collection were less accurate due to the 
distortions discussed in the section on performance analysis. Based on an observation that 
the DOA intersection points of the true DOAs tended to jitter around the true locations, a 
more robust metric which uses the center point of two subsequent DOA intersections was 
developed; this method can provide an estimate of the correct DOA in four collections. 


Conclusion 

It is hoped that the results of this thesis, which are a first step towards deeper investigations 
into the subspace direction-finding algorithms, will enable an expansion of the geolocation 
research which has been on-going at the Naval Postgraduate School to include more anal¬ 
ysis beyond the more traditional two-receive element geolocation methods. The ability to 
accurately utilize MUSIC in moderate to low bandwidth environments, despite it being de¬ 
rived specifically for sinusoidal signals, was also demonstrated in this thesis. Additionally, 
the density of radio frequency (RF) emitters is increasing around the world, and with that 
comes an increase in the likelihood of receiving overlapping signals. Methods of accurately 
resolving and geolocation those overlapping signals are likely to be of interest to the mil¬ 
itary in the future, which means that now is the time for research into subspace direction 
finding and other promising methods of multiple signal geolocation. 



List of References 


[1] Department of Defense Dictionary of Military and Associated Terms, Joint Publication 1-02, 
U.S. Joint Chiefs of Staff, Washington, DC, 8 Nov. 2010, (As amended through 15 Aug. 2014). 

[2] R. Schmidt, “Multiple emitter location and signal parameter estimation,” Antennas and 
Propagation, IEEE Transactions on, vol. 34, no. 3, pp. 276-280, Mar. 1986. 

[3] R. Roy and T. Kailath, “ESPRIT-estimation of signal parameters via rotational invariance 
techniques,” IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 37, no. 7, 
pp. 984-995,July 1989. 

[4] A. Barabell, “Improving the resolution performance of eigenstructure-based direction-finding 
algorithms,” in Acoustics, Speech, and Signal Processing, IEEE International Conference on 
ICASSP ’83., vol. 8, Apr 1983, pp. 336-339. 

[5] B. Friedlander and A. Weiss, “Direction finding for wide-hand signals using an interpolated 
array,” IEEE Transactions on Signal Processing, vol. 41, no. 4, pp. 1618-1634, Apr. 1993. 

[6] M. Doron, E. Doron, and A. Weiss, “Coherent wide-band processing for arbitrary array 
geometry,” IEEE Transactions on Signal Processing, vol. 41, no. 1, pp. 414^17, Jan. 1993. 



THIS PAGE INTENTIONALLY LEET BLANK 


XX 



Acknowledgments 


I would like to start by thanking my advisors, Professor Herschel Loomis and Professor 
Frank Kragh, for their support and direction as I worked to gain an understanding of the 
material and determine where I might, in some small way, be able to add to the knowledge 
base of those who came before me. Your advice has always been exactly what I needed to 
focus my efforts. 

To the faculty, staff, and fellow students of the Space Systems Academic Group, I would 
like to say a sincere thank you for broadening the horizons of this lowly helicopter pilot 
and helping to prepare me for working in the field of space acquisitions. 

Lastly but most importantly I would like to thank my wife for her unceasing love and 
support. You have continually and selflessly picked up responsibilities at home that I have 
had to neglect to focus more and more time on classes, projects, and this thesis. I thank 
God every day for bringing you into my life! 




THIS PAGE INTENTIONALLY LEET BLANK 


xxii 



CHAPTER 1: 
Introduction 


1.1 Background 

With radio-frequency signals, one is frequently only concerned with the signal content, or 
message. There are numerous situations and signal types, however, where descriptive in¬ 
formation about the signal may be more useful. With radar, the timing and doppler return 
of a signal provide information about a target. With beacons, the desire is to determine 
one’s position with respect to one or more fixed emitters, or an emitter’s location with 
respect to a known position. For the military. Joint Publication 1-02 defines signals in¬ 
telligence (SIGINT) as intelligence derived from communications, electronic, and foreign 
instrumentation signals [1]. By examining received radar pulses, determinations can be 
made about the type of radar, its location, and sometimes even the specific piece of equip¬ 
ment; similarly, through evaluation of received communication signals, determinations can 
be made about an emitter’s location and sometimes about the content being transmitted. 
Passive direction finding is the process of estimating an emitter’s direction from the re¬ 
ceiver’s position, and geolocation is then the process of using one or more direction-of- 
arrival (DOA) results in combination with known information about receiver or emitter 
motion and relative geometry in order to determine the emitter’s position in terms of an 
earth-fixed coordinate system. For the military, obtaining location information for emitters 
is a vital part of understanding an adversary’s order of battle and force strength. 

In detecting a signal and obtaining an accurate geolocation, numerous complications can 
arise. Low power and short duration signals are both more difficult to detect. Additionally, 
the task of geolocation is complicated when multiple signals overlap in time, frequency, or 
direction. Airborne and space-based collectors are ideal for conducting SIGINT because 
they provide more area coverage and their motion enables finer location resolution than 
terrestrial systems, but they have limitations on the size and number of receivers which 
can be carried. While the time overlap of short-duration signals are generally unlikely, the 
likelihood increases with increased collector coverage area and in particularly signal-dense 
environments. Many of the current algorithms for separating overlapping signals are based 


1 




on a narrowband assumption or require strict receiver antenna array arrangements, which 
can limit their usefulness [2], [3]. 


1.2 Objectives 

The implementation and performance of the multiple signals classification (MUSIC) sub¬ 
space direction-finding algorithm with a sparse receive antenna array is examined in this 
thesis. MUSIC is based on a narrowband assumption, and it is a goal of this thesis to 
evaluate how wideband the input signals can be before the algorithm breaks down. Basic 
representative types of narrowband and wideband signals are utilized which may overlap 
in time, frequency, and direction and the impact of those overlaps and interactions on the 
ability to accurately determine an emitter’s DOA is evaluated. Specific focus is made on 
the effects of signal frequency, bandwidth, DOA, and signal-to-noise ratio (SNR) on the 
MUSIC output. Additionally, a method of improving the computational performance of 
MUSIC without sacrificing accuracy is presented. Finally, a method of utilizing collec¬ 
tor motion to determine emitter location in the presence of multiple ambiguous DOAs is 
developed for the case of a sparse array with inherent ambiguities. 

This analysis was conducted solely as a software simulation written in MATLAB. In 
addition to implementation of the MUSIC algorithm, test signal generation, a least-squares 
emitter location estimation algorithm, and a signal reconstruction algorithm based on the 
estimated DOAs have also been implemented. The code is included in Appendix B for the 
reader to reference or utilize for future work. 

1.3 Related Work 

The mathematics behind traditional passive direction finding and geolocation have been 
presented in numerous papers and reports. Specifically for this thesis, the clear presentation 
of the methods of geolocation for two receiver elements in reports by Professor Herschel 
Loomis [4] and Michael Price [5] were invaluable to gaining an understanding of the basic 
direction finding and geolocation processes and the uncertainties involved. 

In order to determine the DOA of multiple simultaneous signals, a receiver antenna 
array with three or more receive elements is required. There are two primary methods 
of performing the analysis of such an array—maximum-likelihood and signal subspace 
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methods. The maximum-likelihood estimator is based on a minimization of the likelihood 
function, which is a non-linear optimization problem and is generally computationally in¬ 
tensive to calculate, while subspace methods are based on the eigen-decomposition of the 
signal covariance matrix, for which an estimate is relatively simple to compute [6]. The 
MUSIC method of direction finding developed by Schmidt [2] was one of the first sub¬ 
space methods proposed; since its proposal there have been multiple modifications to the 
MUSIC algorithm to improve resolution or reduce computational complexity, most notably 
Barabell’s root-MUSIC [7]. The advantages of the original MUSIC algorithm are that it is 
relatively simple to understand and can be applied to receiver arrays of arbitrary spacing, 
though it is still computationally intensive because it relies on a search across all possi¬ 
ble signal DOAs. Root-MUSIC avoids this search by examining the roots of the spectrum 
polynomial, but is only applicable for uniformly spaced arrays. Another method proposed 
by Roy and Kailath called estimation of signal parameters via rotational invariance tech¬ 
niques (ESPRIT) [3] is a popular subspace method which is contingent on a matched pair 
of arrays of which the sensors all have identical displacement vectors; like root-MUSIC 
it also avoids a computationally intensive search, but it is a poorer choice for aerospace 
applications due to the requirement for twice as many receivers as the original MUSIC 
algorithm. 

Two drawbacks to all of the subspace methods listed above are that they are non¬ 
coherent, meaning that they will not work if two or more of the received signals are highly 
correlated, and they are designed to work with narrowband signals. The first drawback can 
be overcome via application of spatial smoothing [8], [9]. The second drawback has led 
to the development of various alternative processing methods. Two similar methods which 
derive a composite covariance matrix through a combination of interpolations of the array 
manifold for different frequency bands were developed by Weiss, Doron, and Friedlan- 
der [10], [11]. Any of the subspace methods listed in the previous paragraph can then be 
applied to this modified covariance matrix; these methods also have the added advantage 
of being coherent and able to handle correlated signals. 

In terms of performance analysis, Kaveh and Barabell [12] conducted an in-depth sta¬ 
tistical analysis of the resolution threshold with respect to SNR of two very closely spaced 
narrowband plane waves, and Swindlehurst and Kailath [13] investigated the performance 
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of MUSIC in light of the uncertainties inherent in estimation of the array covariance matrix 
and the array manifold. An analysis of the performance of a sparse array in the presence 
of wideband signals of varying frequencies is not readily apparent in the literature, though 
an understanding of these cases is useful for the implementation of a subspace algorithm 
in an airborne or spaceborne environment and for being able to make informed decisions 
about appropriate selection of a wideband correction method, if one is required at all for 
the target signals-of-interest. 

1.4 Thesis Organization 

Some background information on the importance of emitter location, especially for the 
military, was provided in this chapter. The thesis objective was also presented as was 
related research on which this thesis builds. 

The mathematical algorithms upon which the thesis work relies is discussed in Chapter 
2, Emitter Location Processing. First presented are the basic types of direction finding pro¬ 
cessing with one or two receive antennas, followed by a discussion of geolocation methods. 
The MUSIC algorithm and the underlying concepts of the matrix pseudo-inverse, eigen- 
decompositions, and covariance matrices are then detailed. A subspace visual depiction of 
the MUSIC algorithm is also provided to assist with understanding the algorithm from a 
geometric perspective. 

In Chapter 3, MUSIC Performance Analysis, an in-depth analysis of the MUSIC al¬ 
gorithm’s response to variations in the frequency, bandwidth, DOA, and SNR of a set of 
standardized test signals is provided. In particular, situations where the algorithm breaks 
down and the results are degraded or distorted are pursued and analyzed for their most 
likely causes. 

In Chapter 4, MUSIC Performance Enhancements and Limitations, two enhancements 
to the MUSIC algorithm and some of the limitations inherent in a real-world implementa¬ 
tion are provided. The first enhancement is a method of reducing the computational com¬ 
plexity without loss of resolution, and the second is a method of using collector motion to 
resolve array ambiguities. 

In Chapter 5, Conclusion, the analysis conducted, its significance, and recommendations 
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for future work are discussed. 


In Appendix A, Subspace Processing Specifics, the differences and similarities of time- 
domain and frequency-domain methods of subspace processing as well as considerations 
for the implementation of MUSIC processing at intermediate frequencies are discussed. In 
Appendix B, MATLAB Code, the MATLAB code used to implement the MUSIC algorithm 
and conduct emitter geolocation is provided. 
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CHAPTER 2: 

Emitter Location Processing 


2.1 Basic Direction Finding 

There are four primary methods of passively determining the DOA of signals when only 
one or two receivers are utilized: angle-of-arrival (AOA), interferometry, time difference- 
of-arrival (TDOA), and frequency difference-of-arrival (FDOA). The following paragraphs 
will expand on the processing methodologies for each in two dimensions; the expansion to 
three dimensions is straightforward. The methods are listed below in order of increasing 
processing complexity. 

AOA requires the least amount of signal processing, but it is contingent on the use of a 
high-gain, narrow-beamwidth receiver. The estimated signal DOA is simply the direction 
in which the receiver is pointed, or scanned in the case of a phased array antenna. AOA 
is also unique among the direction-finding methods for only requiring one receiver; it is 
historically referred to as radio direction finding [4]. 

Interferometry could also be called phase difference-of-arrival, because the phase dif¬ 
ference of the signal incident on two or more receiving elements is exploited to determine 
the DOA [5]. A depiction of the two-receive element case in two dimensions, where the 
distance between receive elements (7?i,7?2) in the array is d, is provided in Fig. 1. The 
emitter {E) location is defined in terms of polar coordinates with distance r and angle from 
receiver array broadside 0, which is defined from a common origin point of the array. Ad¬ 
ditionally, the emitter is assumed to be in the far field {r » d) with the signal arriving 
at the receivers as a plane wave. Unless the emitter is directly broadside to the receivers 
(0 = 0°), this plane wave will arrive at one of the receivers before the other, resulting in a 
phase difference (j) at each of the receiving elements. The path difference d which the wave 
travels to reach the second receiver is related to the DOA 0 and the receiver separation d 
by the relationship 5 = Jsin0. From this, the phase difference can be determined as 

^ = 2;rYsin0, (2.1) 

A 


1 




where the wavenumber k = In/X and X is the wavelength; (2.1) also assumes that the 
signal frequeney and wavelength do not ehange during the additional transit time. Deter¬ 
mining the DOA is then simply a matter of solving for 9 in (2.1) [5]. One of the primary 
drawbaeks of this method is that sinee the phase is eyelieal over (j) G [—;r, ;r), the maximum 
unambiguous spaeing between elements over 6 G [—90°,90°] is d = A/2. At separations 
larger than this, the number of ambiguities present will depend ond, X, and 0. 



Figure 1. Two-element interferometer, after Ref. [5]. 


TDOA takes advantage of the same differenee in path length as interferometry exeept 
that the time differenee T is being utilized viee the phase differenee. For TDOA, the path 
differenee is defined as 5 = cT. Solving for 0 direetly as was done in the interferometrie 
ease limits T to a narrow interval in order to avoid phase ambiguities; instead one ean treat 
the emitter’s loeation as an unknown position (xe^ye) and express the path differenee as [4] 


d = 



( 2 . 2 ) 


assuming that the reeeivers are equidistant from the origin and lie on the x-axis. Equation 
(2.2) ean be re-arranged to the form [4] 


2 2 

52 J2_52’ 

T 4 


(2.3) 
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which is recognizable as the formula for a hyperbola. This hyperbola, known as an 
isoehrone, conneets all the possible emitter locations which result in the same TDOA; 
the proper leg of the hyperbola over 0 G (—90°, 90°) ean be determined by the sign of the 
path differenee. An example of a set of isoehrones is presented in Fig. 2. In this manner, 
the only practical limit to receiver spaeing lies in ensuring that both reeeivers ean still re- 
eeive the same signal wavefront at approximately the same power, and reeeivers ean thus be 
spaeed mueh further apart than A/2. The TDOA method of direetion finding is generally 
only applied by itself to signals whieh have eharaeteristies enabling precise time-of-arrival 
determination such as radar pulses. 



Rather than being based on the path differenee, FDOA arises from an artifaet of reeeiver 
motion during signal eollection - Doppler effeet. Since the two reeeivers have different 
veloeities relative to the emitter, a differenee in Doppler frequeney ean be determined and 
exploited to assist in geolocation [5]. Making the assumption that the emitter is not moving, 
we find that the Doppler shift at one reeeiver is the product of the frequency with the ratio 
of the elosing veloeity with the speed of light: A/ = fvdosing/c- The FDOA y is then 
the Doppler difference between the two emitters, y = f{v dosing, i — v dosing, 2 )/c. Sinee the 
closing velocity can be expressed as the dot product of the receiver velocity and the unit 
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vector in the direction of the emitter, in the x-y plane the FDOA can be expressed as [4] 

y-J ( ~ ^ (2 4) 

^ \Vi^e + s)^+y^ ^{Xe-s)^+y^ ) ’ 

under the assumption that the receivers have the same velocity vector and are both on the 
v-axis at distance 5 from the origin, and where and Vy are the closing velocities in the 
X- and y-directions. The curves connecting points of equal FDOA are termed isodops. 
While not obvious from (2.4), this equation defines a leminscate, or “two-leafed rose,” 
when the receivers are moving along the x-axis, which is illustrated by the example set 
of isodops in Fig. 3; when the receivers are moving along the y-axis, a “four-leafed rose” 
shape arises [5]. While a single FDOA isodop does not provide an indication of signal 
DOA like a single isochrone from TDOA or interferometry, it is highly useful in improving 
the position accuracy when used with an isochrone. 

2.2 The Cross-Ambiguity Function 

When either the transmitter or receivers are moving, the Doppler difference in the received 
signals prevents the determination of TDOA via direct correlation, and the cross-ambiguity 
function (CAF) must be used to estimate both TDOA and FDOA simultaneously. The 
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generalized form of the CAP is [14] 


AF(T,r)= r + (2.5) 

Jo 

where 5 i and S 2 are the analytie signals at the two reeeivers, * indieates the eomplex eonju- 
gate, and T and 7 are the possible TDOA and FDOA values over whieh the analysis is eon- 
dueted. The estimated TDOA and FDOA of the eolleeted signal will be where |AF(t, 7 ) | 
has a peak [14]. The TDOA and FDOA results of the CAP improve in resolution with 
inereases in the signal bandwidth and the eolleetion duration, respeetively; thus, the rela¬ 
tively long eolleetion durations possible on eontinuous wave signals also make them prime 
eandidates for analysis via CAP due to the fine FDOA resolution [14]. 

2.3 Position Estimation 

A set of DOA estimates obtained from any of the methods above ean be eombined to 
estimate the emitter’s loeation. In three-dimensional spaee, a minimum of three DOAs is 
required for a position estimate, though two ean be used with an assumption that the emitter 
is loeated on the surfaee of the earth. Use of more DOA measurements will result in an 
over-eonstrained solution; however, due to noise and measurement and timing inaeeuraeies, 
the DOAs will not all interseet at a single point. In sueh a ease a least-squares solution is 
desired, whieh determines the point of minimum residue. Sinee their defining equations 
((2.3), (2.4)) are non-linear, the traditional solution using TDOA, FDOA, or CAP results is 
to derive a linear approximation in the vieinity of an estimated position and use an iterative 
method sueh as Newton-Raphson to eonverge on the loeal minima [4]. 

2.4 Subspace Methods of Direction Finding 

When the possibility exists of q signals being simultaneously ineident on a eolleetor, the 
individual signal DOAs generally eannot be determined unless the eolleetor employs an 
array of p reeeivers, where p > q. The previously diseussed direetion-finding methods 
are designed for use on a single signal ineident on two reeeivers and eannot be direetly 
applied; therefore, alternative methods whieh ean simultaneously proeess information from 
all p reeeivers are desirable. The subspaee direetion-finding methods eomprise one set of 
these alternative teehniques. For a generie situation with q signals and p reeeivers, the 
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sum-of-signals at each receiver can be viewed as the vector multiplication problem 


X = As + w, 


(2.6) 


where x is the p-length vector of received signals, A is the signal steering matrix (often 
referred to as the array manifold), s is the ^-length vector of source signals, and w is the 
p-length additive noise process [15]. For the azimuth-only case, (2.6) can be expanded as 


Xl{t) 


«i(0i) ai(02) ■ 

■■ aiiGq) 


Sl{t) 


w\{t) 


Xlit) 

= 

«2(0l) 



52 (t) 

+ 

wiit) 

(2.7) 

_Xp{t)_ 



^pi^q) 


Sqit)_ 


\Vp{t)_ 



to highlight the structure of the steering matrix, where aj{6i) is the phase shift of the zth 
source signal at the jth receiver. It should be noted that while (2.7) is expressed as a time- 
domain relationship, it can just as easily be expressed in the frequency domain. Addition¬ 
ally, the math is simplified by working with analytic signals; thus, if an arbitrary signal 
Si{t) is a sine wave at frequency ft, it would be represented as Si{t) = 

Since the real portion of all of the signals is utilized, it is simpler to drop the Re{} notation 
and work with just the complex representation, which is known as the analytic signal. In 
this situation, the a(0,) columns of the steering matrix, known as the steering vectors, are 
a p-length vector of the phase shifts applied to each signal relative to the distance of the 
receive elements from a common origin: 


a(0/) = 





^j^ip 



( 2 . 8 ) 


where = kdjsin{6i) for a linear array as defined in (2.1), and dj is the distance from the 
jth receive element to the origin. The phase shift is positive for positive dj and 9 because 
the signal at the jth emitter leads the signal at the origin; in other words the wavefront 
arrives at the jth emitter prior to arriving at the origin. The a{9) steering vectors serve to 
map the signals into p-dimensional space, and the range of x will be limited by the range 
of A, neglecting the effects of the additive noise w. Given that x and A are known and 
s is desired, since p > q the problem is over-determined and a least-squares solution can 
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generally be found, understanding that the result will be noisy sinee w is unknown. Sinee 
6 (and, thus. A) and s are both unknown, a direet least-squares problem eannot be solved. 

Subspaee direetion-finding methods solve this problem by exploiting the eigen-deeomposition 
of the array spatial eovarianee matrix [15]. In general, the eovarianee funetion between 
times ti and t 2 of two signals x{t) and y{t) whieh are modeled as random proeesses is 
defined as [6] 


Cxy itut 2 )=E{ [x{ti ) - E {x{ti )}] [y{t 2 ) -E{y{t 2 )}]*} , (2.9) 

where is the expeetation operator. The eovarianee of the same zero mean signal x{t) 
with itself eollapses to the autoeorrelation funetion. 


Rxxitl,t 2 ) =E{x{ti)x*{t 2 )}. 


( 2 . 10 ) 


Sinee the signals of interest in this direetion-finding problem are sums of sinusoids, 
they are by definition zero-mean. In matrix form over the reeeived signal veetor x then, the 
eovarianee matrix of zero-mean signals ean be defined as 

Rxx=£{xx^}, (2.11) 

where ^ represents the eonjugate transpose, also known as the Hermitian operator. It should 
be noted at this point that C and R are traditionally used by many textbooks to denote the 
eovarianee and eorrelation funetions, respeetively. Most of the subspaee method literature 
appears to use the R notation but defines it as the array eovarianee matrix; sinee zero-mean 
signals are generally assumed, the eovarianee and eorrelation matriees will be equal, but 
the notation ean still be slightly eonfusing. The eonvention utilized by previous subspaee 
papers is used in this thesis. In praetiee, Rxx cannot be determined exactly since the signal 
information is not known for all time, but an estimate can be made over a collected sample 
set as [6] 

Rxx = ^XX", (2.12) 

where the rows of X are the collected samples, making X px N, assuming N samples at 
each receiver. Substituting (2.6) into (2.11), we derive the px p array covariance matrix in 
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terms of the transmitted signal and noise 

Rxx = E {xx^ }='£’{ (As + w) (As + w)^ } = 

£■ {Ass^A^ + w(As)^ + (As)w^+ww^} = , (2.13) 

E { Ass^A^ + ww^ } = ARssA^ + 

where Rss is the source signal covariance matrix and R^^ is the noise covariance matrix. 
In this derivation, the assumption is made that s and w are uncorrelated. If the further 
assumption is made that the noise is additive white gaussian noise (AWGN), then R^w = 
CT^Ip, where aw is the standard deviation of the noise and Ip is the px p identity matrix. 
By construction, Rss is qx q and rank pre- and post- multiplying by A makes ARssA^ 
p X p, but it will still be rank deficient with rank q. The received signal covariance matrix 
Rxx and the AWGN awlp are px p with rank p, and some conclusions can then be drawn 
about the eigen-space of each of the matrices. 

An eigenvalue/eigenvector pair are a scalar and a vector which satisfy the equation 

Av = Av, (2.14) 

where A is an arbitrary square matrix, v is an eigenvector, and A is the eigenvalue. There are 
numerous methods of conducting the eigen-decomposition of matrix which can be found 
in many textbooks and will not be reviewed here. It is often useful to think about eigen¬ 
vectors in a geometric sense: let B be a 2 x A set of (x,y) sampled measurements which 
are randomly distributed over an ellipse, such as a cannon being fired N times and the 
cannonball’s landing position being recorded. The covariance matrix of the data would 
then be the 2 x 2 matrix Rbb This matrix will have two 2x1 eigenvectors which will 
be vectors in the directions of the major and minor axes of the data, and the eigenvalue 
associated with the semi-major eigenvector will be larger than that associated with the 
semi-minor eigenvector. Additionally, a covariance matrix is always Hermitian symmetric 
and positive semi-definite, so the eigenvalues will always be > 0 and the eigenvectors will 
be orthogonal. The eigenvalues will also be equal to the variance of the original data in 
that direction [ 6 ]. If the set of samples in B were all exactly along a line, then one of the 
eigenvalues of Rbb would be zero since the matrix B has no information in the direction of 
the eigenvector associated with the zero eigenvalue (i.e., zero variance); since orthogonal 
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eigenvectors form a multidimensional basis, this means there is no information in one of 
the dimensions. 

Returning to (2.13), since is rank q and is positive semi-definite by construc¬ 

tion, it will have p — q eigenvalues which are exactly equal to zero since a rank deficiency 
indicates that there is no information in one or more dimensions. The AWGN, however, 
notionally exists in all directions equally with p repeated eigenvalues of This means 
that if the eigenvectors of (J^Ip are matched exactly to the eigenvectors of ARssA^, the 
eigenvalues of Rxx will be a combination of the eigenvalues of ARssA^ and G^Ip. Since 
p — q eigenvalues of ARssA^ are zero, p — q eigenvalues of Rxx rnust be equal to be¬ 
cause the only contribution in those dimensions is the AWGN. Additionally, assuming the 
SNR is greater than 0 dB, the q non-zero eigenvalues of ARssA^ are larger than <7^ and 
the contribution from the array covariance matrix will dominate the q largest eigenvalues 
of Rxx- Therefore, the received signal subspace can be thought of as having q signal eigen¬ 
vectors (Es) and p — q noise eigenvectors (E^). This is the basic predicate behind all of 
the subspace direction-finding methods, but they differ in how they exploit this relation¬ 
ship [15]. 

While this derivation was done using AWGN, the same method still applies in the pres¬ 
ence of colored noise. The proof of eigenvalue/eigenvector relationships between Rxx, 
ARssA^, and R^^ in the colored nose case can be shown via the theory of generalized 
eigenvectors [3], and if the noise covariance matrix is known, it can be used to whiten 
the receiver data prior to processing [15] and can be used in estimation of the incident 
signals [2]. 

2.4.1 Subspace Position Estimation 

As discussed in Section 2.3, TDOA and FDOA require an iterative geolocation method such 
as Newton-Raphson because isochrones and isodops are curved contours. The subspace 
methods such as MUSIC are more related to basic interferometry with the direction finding 
output being a DOA. Since a DOA is essentially a vector with unknown magnitude, a direct 
solution can be determined from a set of collections of the signal from the same emitter, 
without having to iterate. 

With only two collections the position estimation is trivial since it is simply the crossing 
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point of two vectors. With three or more collections, the vectors are unlikely to cross at 
the exact same position, and a least-squares estimate must be made. For the more general 
three-dimensional case, the distance of a point from a line can be defined as [16] 


h = 


d X (a —c)| 

m 


(2.15) 


where a is a point on the line, d is a vector along the line, c is the point in question, 
II II indicates the 2-norm, and x indicates the cross product. For this application, a is the 
collector location, d is the DOA vector, and c is the estimated emitter location. A least- 
squares solution for multiple lines can be obtained by solving for the value of c for which 
the derivative of the sum of squares is zero: 


y d jh^) 

h dc 


0 , 


(2.16) 


where n is the total number of collections being used to estimate the position, the notation 
of a derivative with respect to a vector indicates that the derivative should be taken with 
respect to each dimension separately, and 0 is the zero vector. For the two-dimensional 
case, and additionally assuming d is a unit vector, expanding (2.15) results in 


h — Cxdy Cydx Clxdy Hydx- 


(2.17) 


By squaring (2.17), taking the derivative, and then summing the results for all n collections 
as specified by (2.16), we can convert the problem into an AX = B type matrix equation 
where c can be solved for directly: 


n 


I 




n 


I 



dxidyi 


dxidyi 




(2.18) 


2.4.2 Signal Reconstruction 

Once an estimate has been made of signal DOAs, the original signals can be reconstructed 
by using (2.8) and (2.1) to recreate the steering matrix and then solving for s in (2.6) as 
s = A^^x. Generally, the steering matrix will not be square, however, so its inverse cannot 
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be taken directly. Instead, the pseudo-inverse can be calculated. For an mx n matrix A 
where m> n and the columns of A are linearly independent, the Moore-Penrose pseudo¬ 
inverse is defined as [17] 

\+=[A^Ay^A^. (2.19) 

Since p > q and the steering matrix spans the signal subspace, its columns are linearly in¬ 
dependent and this definition of the pseudo-inverse applies. Therefore, the original signals 
can be estimated as 

s = A+x. (2.20) 


2.5 MUSIC 


The MUSIC algorithm was proposed by Ralph Schmidt in 1979 [2] and was one of the 
earliest subspace methods. MUSIC develops from a recognition that since the eigenvec¬ 
tors of the array covariance matrix are a basis and are orthogonal to each other, the noise 
eigenvectors will be orthogonal to the signal eigenvectors. Since the steering vectors must 
by definition lie in the signal subspace, they will also be orthogonal to the noise eigenvec¬ 
tors. Recall from (2.8) that, given a known array center frequency and layout, the steering 
vector a is only a function of the signal DOA 0. The MUSIC algorithm is then a search 
over all 9 in the visible range for the steering vectors which are most orthogonal to the 
noise eigenvectors of the array covariance matrix, resulting in peaks in the MUSIC DOA 
spectrum P [2]: 


PMUSicid) 


a^(0)a(0) 

a"(0)EwE^^a(0)' 


( 2 . 21 ) 


Note that the normalization in the numerator is only required when the steering vectors are 
not already normalized to magnitude 1. An example of the MUSIC spectrum using the 
parameters from the following section is plotted in Fig. 4. Note that the y-axis is labeled 
“Orthogonality”; since the denominator of the MUSIC algorithm is effectively an inner 
product, the MUSIC response at each DOA is an indication of how orthogonal the steering 
vector at each DOA is to the noise eigenvector(s). At exactly the correct DOA, assuming 
no receiver noise, the MUSIC response should be infinite because the inner product will be 
zero. Since the plots can have such a large dynamic range (theoretically infinite), the y-axis 
is plotted in decibels. 
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Figure 4. Sample MUSIC DOA spectrum for emitters at 74° and 116°. 

In any real-world implementation with a finite number of samples with which to esti¬ 
mate the array covariance matrix, Oy, will not be exactly equal in all directions, and the 
resulting noise eigenvectors will not be exactly orthogonal to the original signal subspace. 
This will lead to slight errors in the MUSIC spectrum peaks for any particular collection 
period. MUSIC is asymptotically unbiased, however, so the spectra averaged over multiple 
collections will converge to the exact solution [2]. 

2.5.1 MUSIC Example 

For the purposes of this example, let the receive elements be isotropic and arranged as the 
points on an equilateral triangle with receiver spacing of A/2 and the origin at the center, 
as indicated in Fig. 5, where A is based on the array center frequency of 200 MHz. Emitter 
1 is transmitting a continuous cosine wave at 200 MHz and 6i = 116°, and Emitter 2 is 
transmitting a 205 MHz cosine wave with 02 = 74°, and we assume there is no additive 
noise in the receiver. The collection period is 100 ns and the sampling frequency is 800 
MHz. 

Since this array layout is not linear, the receiver array does not have a broadside, so 
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the ^ij terms in (2.8) must be defined slightly differently. For this arrangement let = 
kdj cos aij, where dj is the distanee to the reeeiver from the origin and atj is the angle 
between the jth reeeiver, the origin, and the ith emitter, whieh is also equivalent to the 
differenee between the angle to the emitter and the angle to the reeeiver relative to the x- 
axis. All three reeeivers are X/{2^/3) from the origin, Ri is at 0° from the v-axis, R 2 is 
at 120 ° from the x-axis, and R^ is at — 120 ° from the x-axis, and the signal arrival angles 
relative to the x-axis have already been speeified. The produet 


^ ^ 2;r A 7t 

and the relative angles will equal 



' 14--0° 

1 

0 

0 

1 


■ 740 

116°' 

a = 

74-120° 

116-120° 


-46° 

-4° 


74+120° 

116+120° 


194° 

236° 


( 2 . 22 ) 


(2.23) 
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The steering matrix is then 



- ^7^cos74° 

^y^cosll6°- 


■ ^JO.500 

g-70.795- 

A = 



= 

1.260 

1.809 


;■+= cos 194° 

^y^cos236° 


g-;i.760 

g-7l.014 


(2.24) 


Since w = 0 and s = 


J(27r/lf) pj{27tf2t) 


iT 


, where /i = 200 MHz and /2 = 205 MHz, the 


received signal vector can now be defined: 


X = As + w = 


+ 0 . 500 ) _|_ ^j2Kf2t-Q.195) ' 
^j{27tht+\.26Q) j^^j{2nf2t+l.m) 
^j{27tfit- 1 . 760 ) ^ gj{2Kf2t- 1 . 014 ) 


(2.25) 


The sample estimated array covariance matrix Rxx must be determined numerically; 
using MATLAB to generate the 3 x 100 received signal matrix and 3x3 covariance matrix 
results in 


Rxx 


0.7753 

-0.0901 + jO.7374 
0.2370-jO.6537 


-0.0901-jO.7374 
2.6818 

-2.7092-jO.6237 


0.2370 + jO.6537 
-2.7092 + jO.6237 
2.8917 


(2.26) 


which, after an eigen-decomposition, has the eigenvalues and eigenvectors which are listed 
in Table 1. Note that the third eigenvalue is zero because there is no signal (or noise) present 
in this dimension; therefore, the third eigenvector is the only noise eigenvector. 

Table 1. Circular array eigenvectors and eigenvalues. 


eigenvalue 

5.7753 

0.5734 

0 


0.0679 + 70.1873 

-0.3326-70.9175 

0.0304 + 70.0838 

eigenvector 

-0.6637 + 70.1528 

-0.0706 + 70.0163 

0.7101-70.1635 


0.7046 

0.2059 

0.6791 


The MUSIC DOA spectrum can then be obtained using (2.21); the result has already 
been presented in Fig. 4, and the peaks are at the source signal DOAs of 74° and 116°. 
The additional hump at 260° is a near orthogonality which occurs at the mirrored DOA 
angles. Because the plot peaks approach infinity but will not typically equal it due to 
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additive noise, limitations in the angular step size whieh is searehed, and rounding errors 
during proeessing, the DOA speetrum is typieally plotted in deeibels so that the various 
peak heights are visually similar. 

The sharp-eyed reader may have notieed a diserepaney in the above example, in that 
the wavelength used to eonstruet the steering veetor matrix in (2.24) is the array eenter 
wavelength of 200 MHz for both signals, despite the faet that one of the signals is not at 
the eenter frequeney. This arises because the exact source signal frequencies are not known 
a priori to conducting the subspace analysis to separate the signals, and one must assume 
the array center frequency. The resulting DOA determined for the non-center frequency 
signal will have an induced error which is proportional to the frequency offset. Methods of 
correcting for this offset and for the effects of wideband signals are the focus of the rest of 
this thesis. 

2.5.2 MUSIC Visualization 

As alluded to in the eigen-space discussion, the relationship between x. A, and s can also be 
viewed geometrically. If the set of all possible steering vectors (0=0: 360°) is determined 
and plotted in p-space, it will transcribe a continuous ^-dimensional object. For p = 3 and 
q = 2 this will be a line in three dimensional space; for p = 4 and q = 3 this will be a 
surface in four dimensional space. If the array is linearly arranged, the steering vector 
matrix will transcribe two overlapping 180° arcs; if it is not linearly arranged (i.e. two- or 
three-dimensional), it will transcribe a closed loop. 

If the signals are adjusted to be real vice analytic, the example above can be plotted in 
three-dimensional space, where each dimension represents one of the receiver elements. 
By doing this, the signal and noise eigenvectors can be visually distinguished, and the 
orthogonality of the steering vectors in the signal DOAs to the noise eigenvector can also 
be visually ascertained. 

The results are plotted in Figs. 6 and 7. Even though the source signals are real, since 
the steering vectors are complex phase terms, the received signals are complex, as are the 
covariance matrix and eigenvectors, and the problem has to be separated into real and imag¬ 
inary portions to be depicted . In both of these figures, the gray portions are the received 
signals x, and the blue lines are the eigenvectors of Rxx- From visual inspection, two of 
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Receiver 3 




Receiver 1 Receiver 2 

Figure 7. MUSIC geometric visualization - imaginary portion. 
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the eigenvectors in both the real and imaginary portion plots lie in the signal subspace and 
the third is normal to the signal subspace; this third eigenvector is the noise eigenvector. 

The black loop in both figures is the continuum of all steering vectors, and the red lines 
indicate the steering vectors a(0i) and a ( 62 ) at the signal DOAs. By visual inspection 
the red steering vectors lie in the gray signal subspace and are orthogonal to the noise 
eigenvector. 

What the relationship looks like for the real portion of the signals when noise is in¬ 
cluded at a 10 dB SNR is shown in Fig. 8. Because of the additive noise, the received 
signal subspace is no longer only a two-dimensional plane. However, there are two domi¬ 
nant dimensions which equate to the two largest eigenvalues, and the subspace dimension 
with the smallest eigenvalue is still the noise subspace. The orthogonal relationship of the 
steering vectors to the noise eigenvector remains despite the noise. 



Figure 8. MUSIC geometric visualization - real portion including noise. 


23 





2.6 Emitter Location Processing Review 

In this chapter, the principles behind various direction finding and geolocation methods 
were presented, with particular emphasis on explaining the derivation and providing an 
example of the MUSIC subspace algorithm. In the next chapter, the MUSIC algorithm is 
investigated to determine its performance characteristics for a sparse receiver array in the 
presence of signals with varying frequency, bandwidth, DOA, and SNR. 
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CHAPTER 3: 

MUSIC Performance Analysis 


3.1 Analysis Setup 

For the purposes of this analysis, a linear array of six reeeivers on the ;c-axis with a visual 
range of —90° to 90° relative to array Foresight is modeled. Additionally, it is assumed that 
all of the reeeive elements are isotropie. The array eenter frequeney fc is 1 GHz, the seareh 
bandwidth is 500 MHz to 1500 MHz, and sampling is eondueted at 4 GHz. It should be 
noted that this analysis is eondueted at radio frequeney (RF); see Appendix Seetion A.3 for 
a diseussion of implementing MUSIC at intermediate frequeney (IF). This reeeiver array 
is also referred to as the eolleetor. 


3.1.1 Array Configuration 

In order to ensure low sidelobes and no ambiguities, the ideal reeeiver spaeing would be a 
uniform spaeing of one half wavelength at the desired eolleetion eenter frequeney. With six 
reeeive elements, sueh an array has a very short baseline of only two and a half wavelengths 
and a relatively wide beamwidth. The direetion finding resolution of an array varies with 
the proeessing method and the signal parameters, but in general it is proportional to the 
beamwidth and inversely proportional to the baseline; therefore, a large baseline is usually 
desirable. With uniformly spaeed elements, however, a large baseline means that many 
reeeivers are required, whieh inereases the weight and design eomplexity. Additionally, 
reeeive element size may prevent the plaeement of elements as elose as A/2. Airborne and 
spaeeborne applieations tend to be weight and spaee limited, so for those applieations a 
sparse layout may be more desirable. The tradeoff with a sparse array is that ambiguities 
may exist and the sidelobes are higher than with a uniformly spaeed array. These trade¬ 
offs ean be mitigated somewhat through eareful design, sueh as by aeeepting ambiguities 
outside of a partieular limited visible range. 

The sparse layout used throughout the rest of this ehapter is 


d 


0 0.5 1.5 3.5 6 10 


wavelengths, 


(3.1) 
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where J,- represents the loeation of a reeeiver relative to the origin; defining the first reeeiver 
at the origin serves to help simplify the ealeulations. For the same number of elements, this 
array has a four times increase in baseline over the uniformly spaced array and a narrower 
main beamwidth. The sidelobes for this sparse array are significantly higher than for the 
uniformly spaced array, however, as indicated in Fig. 9. In fact with this spacing scheme, 
there are no ambiguous lobes, but all sidelobes are only ~6 to 8 dB down from the main 
beam. A depiction of the resolution difference between the uniform and sparse arrays using 
the MUSIC algorithm is shown in Fig. 10. In this example there are two sinusoidal signals 
at a 2° separation relative to the collector, but the uniform array (Fig. 10(a)) cannot resolve 
them both and its MUSIC spectrum only has one peak, while the MUSIC spectrum of the 
sparse array (Fig. 10(b)) clearly has two peaks. The reason is that the larger baseline of the 
sparse array results in a narrower MUSIC response. Additionally, it is important to point 
out that what appears as noise in Fig. 10(b) is actually the response of the array’s sidelines 
to the signals present. 



Figure 9. Array patterns with six isotropic elements. 

At this point it should also be noted that in general, results are plotted in sinespace vice 
degrees, where sinespace is defined as the sine of the angles in the visible range. This is 
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(a) Uniform array, A/2 spacing. (b) Sparse array of (3.1) 

Figure 10. Six-element array resolution for 2° emitter separation. 


because the MUSIC response of a signal is the same width across all of sinespace, whereas 
the MUSIC responses plotted in degrees widen out as the signals are further from boresight. 
Plotting in sinespace ensures that the primary distortions will be from signal bandwidth and 
array sparsity. A sinespace step size of 0.001 is used for all of the scenarios in this chapter. 

3.1.2 Emitter Parameters 

The incident signals being modeled are variations of three basic signal types: pure sinusoid, 
linear frequency modulation (LFM), and quadrature phase-shift keying (QPSK). All LFM 
signals are modeled as an “up-chirp” centered on a desired frequency, and all QPSK sig¬ 
nals utilize rectangular pulse shapes with phase terms generated using the randi function 
in MATLAB. It should be noted at this point that pulsed QPSK modulation is not utilized in 
phase-coded radar systems. This signal type has been chosen for this thesis because it has 
a bandwidth profile which provides a reasonable approximation to a bi-phase or polyphase 
coded radar pulse as well as a short segment of a phase modulated communications sig¬ 
nal. Radar terminology for a phase-coded signal will be used when describing the QPSK 
signals: a code is the full set of phases transmitted during the QPSK pulse, and a sub¬ 
code is an individual phase term transmitted during one subpulse [18]; in communications 
terminology the subcode would be a channel symbol [19]. 

Emitters are positioned across the visible range at distances between 100 and 200 nmi. 
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and are assumed to be stationary. All of the emitted signals are modeled as pulses which 
arrive at the array at time zero and last for 100 ns. The emitter parameters being varied in 
this analysis are frequency, bandwidth, DOA, and SNR. For the purposes of simplifying 
signal creation and processing, all signals are generated as analytic signals in the time 
domain. 

3.1.3 Collector Motion 

When multiple collections are simulated, the collector is modeled as being mounted on the 
side of an aircraft which is moving in the positive x direction at 300 knots. The altitude 
is not being modeled as this analysis is conducted in only two dimensions; therefore, the 
emitters are assumed to be in the same plane as the collector. For an aircraft at altitude this 
could be the slant plane, though the intersection of that plane with the earth is not being 
utilized as a constraint in the position estimation which was delineated in Section 2.4.1. 

Collections are made every two minutes, and all signals being modeled are assumed to 
be received at every collection time. Doppler shift due to collector motion is not included 
in the model, and the collector is assumed to be stationary for the duration of a collection. 
This first assumption simplifies the association of signals over multiple collections, and the 
latter two assumptions are reasonable based on the short collection durations. 

3.1.4 Processing Method and Results Presentation 

For this analysis the frequency domain form of (2.6) was utilized (see Appendix A.2). 
Working in the frequency domain provides some advantages because the signal frequencies 
are more readily accessible without additional processing and band-limiting can be utilized 
to reduce computation and noise influence. The source signals are generated in the time 
domain, and the discrete Fourier transform (DFT) (specifically MATLAB’s implementation 
of the fast Fourier transform (FFT)) is then used to convert the received signals to the 
frequency domain for all further processing. 

An overhead “map view” of the collection plane may also be included with the results, 
especially when emitter motion is modeled. In these views the collector position(s) are 
indicated by a black circle and the actual emitter positions by a black cross. Estimated 
DOA lines colored to correspond to each received signal are plotted extending from the 
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collector positions and, when applicable, black diamonds are used to indicate the least 
squares estimated emitter positions. 

3.2 Frequency Effects 

In order to gain an understanding of the effects of frequency on the subspace DOA spec¬ 
trum, three sinusoidal signals at the frequencies and locations in Table 2 were modeled 
over five collections. No receiver noise was included, ensuring that only the frequency 
variations would affect the results. 

Table 2. Signals utilized in modeling frequency effects. 


Signal Number 

Frequency (MHz) 

location (nmi) 

1 

1000 

[-50, 90] 

2 

1323 

[-40, 190] 

3 

811 

[-100, 110] 


The results are presented in Fig. 11. The MUSIC results for all five collections are 
overlaid in Fig. 11(c) so that the relative change in emitter locations can be seen. From 
Fig. 11(d), it can be seen that the DOA lines only intersect in the close vicinity of one of 
the three actual emitter locations - Emitter #1, the one at the array center frequency. 


As initially mentioned in Section 2.5.1, this direction finding and geolocation error 
arises because the center frequency was assumed when constructing the steering vectors 
for use in (2.21). There exists, however, a relatively simple solution which requires under¬ 
standing of the origin of the steering vectors. Without loss of generality, the solution can be 
found by examining the noiseless relationship of the fth emitter with the yth collector from 
( 2 . 8 ): 


Xij = Sie 


j27t-^ sin6,' 


(3.2) 


The effect of the exponential term is simply a phase shift, and the received signal magnitude 
Si will be the same regardless of the DOA. For a particular collection, Xij, su and the 
magnitude of the phase shift are fixed. Note that there are essentially two unknown 
variables in the phase shift, however: the wavelength A,- and the DOA 0,-, which means that 
there is a dependent relationship between the wavelength, the DOA, and the phase shift and 
that there are an infinite number of wavelength/DOA pairs which taken together provide the 
correct phase shift. In order to solve for the DOA using MUSIC, an assumption was made 
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(a) Transmitted signals; first 40 ns. 



(c) MUSIC DOA spectrum. 

Figure 11. Results of MUSIC processing on 
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(b) Received signals; first 40 ns. 



(d) Map view of DOAs. 
narrowband signals with different frequencies. 


about the wavelength; therefore, if the wavelength is ineorrect then so is the resulting DOA. 
The value of the phase shift is fixed, so an equality ean be configured to solve for the actual 
DOA as 

dj dj 

(pij = 2 k Y sin Oi^MUSic = -sin 0,-ev, (3.3) 

4-c ^i,est 

which can be further reduced to 


fc sin Oi^MUSIC = fi,est sin Oldest: 


(3.4) 
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where Xc is the wavelength at the array’s eenter frequeney at whieh the MUSIC algo¬ 
rithm was run, Xi^est and fi^est are the zth signal’s estimated wavelength and frequeney, and 
QiMUSic is the DOA of the zth signal as determined by the MUSIC algorithm. There are 
still two unknowns: finest and digest- Sinee the signal parameters of Si ean be determined via 
signal reeonstruetion, (Seetion 2.4.2), the eenter frequeney of the signal ean be estimated 
and the true DOA ean be estimated as 


6i,est = aresin 


fc . a 

-— sm Oi^MUSic 

Ji,est J 


\ 


(3.5) 


Attempting to aeeurately reeonstruet the signal with an ineorreet DOA may seem bound 
to result in errors; however, it ean easily be proved to provide an aeeurate result. Reeall 
that signal reeonstruetion involves an estimation of the steering veetors for eaeh signal, and 
that the steering veetors are simply the eolleetion of phase shifts applied to eaeh signal at 
eaeh reeeiver. As stated previously, the magnitude of the phase shifts are already known, at 
least to the aeeuraey that the DOAs ean be extraeted from the MUSIC DOA speetrum. 

It is important to note that the signal frequeney ean only be estimated from the reeon- 
strueted signal due both to the nature of the DFT and the presenee of noise in an aetual 
eolleetion. A method for estimated the average frequeney in a signal’s diserete frequeney 
response follows. First, normalize the magnitude of the reeonstrueted signal’s frequeney 
speetrum and filter to the desired seareh bandwidth. Next, reduee the impaet of noise and 
the signal not being eentered in the seareh bandwidth by setting the lowest 10% of the fre¬ 
queney speetrum values to zero. Finally, a weighted mean frequeney ean be determined 
by taking the inner produet of the veetor of DFT frequeney bins with the normalized and 
truneated frequeney speetrum and then dividing by the sum of the normalized and truneated 
frequeney speetrum. The reeonstrueted frequeney speetrum for the first eolleetion is pro¬ 
vided in Fig. 12(a), and estimated frequeneies for eaeh of the five eolleetions are presented 
in Table 3. By eomparison with the original frequeneies listed in Table 2 it ean be seen that 
this method works quite well for these narrowband signals. 

Using (3.5) with the values in Table 3, we ean ealeulate frequeney-eorreeted DOAs 
for eaeh eolleetion. A map view with the original DOAs as dashed lines and the eorreeted 
DOAs as solid lines is illustrated in Fig. 12(b). Additionally, the estimated emitter positions 
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(a) Reconstructed signal frequency spectrum, (b) Map view with frequency-corrected DOAs 

and estimated emitter locations. 

Figure 12. MUSIC DOA frequency correction. 


Table 3. Estimated signal frequencies. 


Collection 

/i (MHz) 

fi (MHz) 

h (MHz) 

1 

1000.0 

1323.2 

810.1 

2 

1000.0 

1323.0 

810.1 

3 

1000.0 

1322.2 

810.2 

4 

1000.0 

1322.2 

810.2 

5 

1000.0 

1323.2 

810.0 


have been plotted (diamonds), and one can see that they align very well with the actual 
positions (crosses). 


An interesting artifact of this frequency effect is that MUSIC can easily resolve two sig¬ 
nals which emanate from the exact same location as long as they are different enough 
in frequency. The MUSIC DOA spectra for two QPSK signals at frequencies of 900 
and 1100 MHz which are co-located at —10° from collector boresight are illustrated in 
Fig. 13(a). Both signals also have a bandwidth of 20% fc and a 10 dB SNR. There are 
clearly two distinct peaks in the MUSIC spectra for each collection, but from the map view 
(Fig. 13(b)) it can be seen that the corrected DOAs overlap with estimated locations very 
close to the correct position. It should be noted that this effect also depends on the angle 
from array boresight; if the signals are at boresight, there will be zero relative phase shifts 
between the elements, and the frequency difference will have no effect and there will be a 
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single MUSIC peak. 



(a) MUSIC DOA spectra. (b) Map view. 

Figure 13. Two 10 dB SNR QPSK signals with center frequencies of 900 MHz and 1100 MHz 
and the same source location. 


This frequeney/DOA dependence can also have the alternate effect of causing two sig¬ 
nals with otherwise sufficient angular separation to overlap because of their differences in 
frequency, making them unresolvable in that collection. 

3.3 Bandwidth Effects 

For the initial bandwidth analysis, all signals were centered at the array center frequency 
of 1 GHz, and bandwidths of 50 MHz, 100 MHz, 200 MHz, 300 MHz, and 400 MHz were 
tested, which are respectively 5%, 10%, 20%, 30%, and 40% of the center frequency fc. 
Additionally, all signals were set at a DOA of 15° with no receiver noise. Only the LFM and 
QPSK signal types were implemented since the bandwidth of the sinusoidal pulse is fixed at 
Bp= I /tp, where tp is the pulse duration. For this analysis, tp = 100 ns and Bp = 10 MHz, 
which is approximately equal to the width of one FFT bin. 

The bandwidth of an LFM signal is being defined as = fhi — fio^ where //„ is the 
highest instantaneous frequency transmitted and fio is the lowest instantaneous frequency 
transmitted. The bandwidth of a phase coded signal such as QPSK is being defined as 
Bqpsk = l/7j, where is the subcode duration [19]. Both of these are noise-equivalent 
bandwidth definitions assuming an infinitely long pulse; the true bandwidths as viewed on 
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a spectrum analyzer would be slightly larger due to the pulsed nature of the signals being 
modeled. For a pulse duration of 100 ns, however, this additional bandwidth is a negligible 
fraction of the carrier frequency as noted in the previous paragraph. 

The MUSIC DOA spectra for a single LFM and QPSK signal at different bandwidths 
are presented in Fig. 14. The change in MUSIC response for both signal types appears 
very similar with a relatively linear reduction in peak amplitude inversely proportional to 
the change in bandwidth. 



Angle off boresight (sinespace) Angle off boresight (sinespace) 


(a) LFM signal. (b) QPSK signal. 

Figure 14. Effect of bandwidth on MUSIC spectrum. 

This is likely because even though the bandwidths are the same, a phase coded signal 
has a sine-shaped frequency spectrum where most of the energy is concentrated near the 
center frequency, whereas the energy in a LFM signal is evenly spread across the bandwidth 
since the frequency spectrum is approximately rectangular. 

The effect of bandwidth still appears relatively minor, however, since a peak still re¬ 
mains at the proper DOA even if it is reduced in amplitude. Once multiple received signals 
are evaluated, signal bandwidth can have more noticeable effects on the MUSIC spectrum. 

One might believe that signals with overlapping bandwidth would lead to problems with 
determining an accurate DOA, but that appears to not be the case in at least a basic two- 
signal scenario. The MUSIC DOA spectrum and the transmitted signal frequency spectrum 
for two QPSK signals at —24°/100 MHz bandwidth and 15°/150 MHz bandwidth, both at 
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the array center frequency and with 10 dB SNR, are illustrated in Fig. 15. Despite the 
complete bandwidth overlap, the signals are both easily resolvable and reconstruction is 
fairly good, especially close to the signal center frequency. 







(a) MUSIC DOA spectrum. (b) Frequency spectra. 

Figure 15. Effect of signal overlap on MUSIC DOA resolution and signal reconstruction. 


It has already been shown that the sparse array being utilized can resolve narrowband 
signals with a small angular separation. If those signals are wideband, MUSIC can still 
resolve them, but other issues can arise even with a large frequency separation. Two ex¬ 
amples utilizing a 750 MHz LFM at —30° and a 1000 MHz QPSK signal at —24° with no 
additive noise are provided in Fig. 16. Note that with the frequency adjustment in (3.5), the 
MUSIC algorithm “sees” the LFM signal at —22°, so this scenario has an effective signal 
separation of 2°. In Fig. 16(a) the signals have a bandwidth of 5% fc, and in Fig. 16(b) 
the signal bandwidths are 10% fc- From an examination of both figures, it can be seen 
that the signals have been blended in reconstruction, and the effect worsens as the relative 
bandwidth increases. This blending results from the various signal frequency bins being 
DOA shifted different amounts relative to the signal center frequency; since the signals are 
so close in angular separation, there is a crossing of frequency bins between signals which 
cannot be perfectly untangled using (2.20). 

If a third signal with particular parameters is added to the previous 10% bandwidth 
scenario, another distortion occurs. One might expect that adding a signal with a large 
angular separation would have no impact and would be easily resolvable, but that is not the 
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(a) MUSIC DOA spectrum. (b) Frequency spectra. 

Figure 16. Transmitted and reconstructed frequency spectra of wideband signals with 2° angular 
separation. 


case if the new signal overlaps in frequeney with one of the prior two. For this seenario a 
1030 MHz LFM signal at 32° has been added, with varying bandwidth of 50 and 200 MHz 
(5% and 20% fc). The results are presented in Table 4 and Fig. 17. When the bandwidth 
of the third signal is inereased to the point where it eompletely overlaps the QPSK signal, 
the MUSIC response beeomes distorted. The two elosely spaeed signals are now no longer 
resolvable as separate signals. Additionally, the MUSIC response to the third signal is no 
longer the sharp peak it was with the partial overlap in bandwidth. From Table 4 it ean be 
seen that the DOA estimations from the MUSIC speetrum in the fully overlapped ease have 
been shifted from their eloser-to-true estimations in the partially overlapped ease. 

Table 4. Effect of signal bandwidth on DOA with multiple closely spaced and frequency over¬ 
lapped signals. 


Signal 

1000 MHz QPSK 

750 MHz LFM 

1030 MHz LFM 

True DOA 
(frequeney adjusted) 

-24.0° 

-22.0° 

33.8° 

MUSIC DOA with 
partial overlap 

-24.5° 

-22.1° 

33.8° 

MUSIC DOA with 

-22.6° 

32.9° 

full overlap 
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(a) MUSIC DOA spectra. 


(b) Frequency spectra of transmitted signals. 


Figure 17. Effect of signal bandwidth on closely spaced and frequency-overlapped signals. 


3.4 DOA Effects 

It has already been noted that resolution deereases with inereased DOA beeause the phase 
differenee is a funetion of the sine of the DOA; as the DOA approaehes 90° the DOA is 
ehanging mueh faster than the phase. For example, a sinespaee step size of 0.01 at array 
broadside represents an angular step of 0.57°, while at 90° off boresight that same sinespaee 
step represents an angular step of 8.11°. The DOA ean have other affeets on the MUSIC 
response as well. For this analysis, signals were generated at 1, 15, 30, 45, and 60° with no 
noise. Signals are eentered at the array eenter frequeney, and LFM and QPSK signals have 
a bandwidth of 20%. A depietion of the responses for a sinusoid is provided in Fig. 18. The 
peak amplitudes vary slightly, but that is only due to the diseretization of the angular step 
impaeting the exaetness of the orthogonality response as eertain steps are slightly eloser to 
the aetual DOA than others; reeall that the ideal peak amplitude is infinite. 

The MUSIC response for a LFM signal is provided in Fig. 19(a). For this signal type 
the response amplitude degrades fairly signifieantly, and the response widens as the angle 
from array boresight inereases. The MUSIC response for a QPSK signal is illustrated in 
Fig. 19(b). Here, the QPSK MUSIC peak amplitudes degrade more slowly than those of a 
LFM signal, unlike in the bandwidth analysis above. This differenee is likely due to the faet 
that a phase-eoded signal has more energy eoneentrated near the array eenter frequeney. 
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Figure 18. Effect of arrival angle on MUSIC spectrum with sinusoidal signal. 


while a LFM signal has its energy spread more evenly over the entire bandwidth, so the 
additional distortion that oeeurs with inereased DOA has less effeet on phase-eoded signals 
sueh as QPSK. 



(a) LFM signal. 


(b) QPSK signal. 


Figure 19. Effect of DOA on MUSIC spectrum with wideband signals. 


That a degradation occurs at all is related to the frequency effects analysis of Section 
3.2. As the DOA increases, the relative additional path distance d between elements also 
increases, which results in an increased phase difference (j) for frequencies not at the signal’s 
center frequency. It is useful to think of each FFT frequency bin as affecting the MUSIC 
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response independently; bins above the signal eenter frequeney map to higher DOAs and 
bins below the signal eenter frequeney map to lower DOAs. The response spreads out over 
a larger set of DOAs and also deereases in amplitude beeause there is less energy in eaeh 
bin as bandwidth inereases. 


3.5 SNR Effects 

For this portion of the analysis, SNRs of 100, 10, 3, 0, —3, and —10 dB were utilized with 
Monte Carlo simulations of 100 trials. Signals are eentered at the array eenter frequeney 
with a DOA of 15°. The average MUSIC response at varying SNR for a sinusoidal signal 
with the six-element sparse array of (3.1) is illustrated in Fig. 20. A degradation of the peak 
amplitude with inereased noise ean be seen, similar to the effeet due to bandwidth. In this 
ease the likely reason is that sinee additive noise serves to rotate the signal eigenveetors 
of the eovarianee matrix away from their true loeations, on average the steering veetors 
will be somewhat less orthogonal to the noise eigenveetors as the SNR deereases. Only the 
sinusoidal signal results are presented beeause the degradation rate is very similar for all 
signal types. 
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Figure 20. Effect of SNR on MUSIC spectrum with sinusoidal signal; average of Monte Carlo 
simulation with 100 trials. 
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The standard deviation Odoa and mean niooA of the estimated MUSIC peak location 
across the 100 trials are listed in Table 5; as the SNR decreases, the standard deviation 
increases, which should be expected. The mean DOA estimations remain very close to 
the true DOA of 15°, however, confirming that MUSIC is indeed an unbiased estimator. 
It should be noted that even the mean of the nearly noiseless case of 100 dB SNR is not 
exactly at the correct DOA of 15°; this is due to the finite sinespace step size of 0.001. In 
the vicinity of the correct DOA, the two closest angular values are 14.95° and 15.01°; thus, 
it can been seen that the algorithm converges on the closest angular step and not necessarily 
the exact DOA. This effect will add to the geolocation uncertainty but should average out 
over multiple collections at different DOAs. 

Table 5. DOA statistics for varying SNR. 


SNR (dB) 

100 

10 

3 

0 

-3 

-10 

(^DOA 

1.07 X 10^^4 

0.015 

0.033 

0.050 

0.080 

0.234 

t^DOA 

15.01 

15.01 

15.00 

15.01 

15.01 

15.02 


The role of the number of array receiver elements was also investigated to see if more 
receivers would result in an averaging effect which would improve the MUSIC peak am¬ 
plitude for low power signals. A Monte Carlo simulation with 100 trials was run for a 
uniformly spaced array with 2, 3, 4, 5, and 6 elements and a sinusoidal signal at —10 dB 
SNR, with the averaged results presented in Fig. 21. The baseline was held constant at 2.5 
wavelengths to ensure that changes in baseline did not affect the results. This does produce 
ambiguities at the lower receiver numbers, but the ambiguous results are outside of the 
0-0.5 sinespace range which is plotted and can be ignored for this scenario. The relatively 
constant peak amplitude is somewhat surprising, as one might expect that by averaging 
across more receivers the effective SNR would improve. 

The standard deviations and means at each array size are listed in Table 6. The means 
vary somewhat randomly, indicating that the estimated DOA accuracy is more a function of 
the baseline, which is constant in this scenario. The standard deviations steadily decrease 
with array size, indicating that the precision over multiple trials does improve with the 
number of receivers. 
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Figure 21. Effect of the number of array elements on the MUSIC spectrum with a sinusoidal 
signal at —10 dB SNR; average of Monte Carlo simulation with 100 trials. 


Table 6. DOA statistics for varying array size. 


Array size 

2 

3 

4 

5 

6 

Odoa (deg) 

1.69 

I.3I 

1.22 

1.02 

0.85 

mooA (deg) 

15.12 

15.15 

15.01 

14.88 

15.12 


While an averaging effect across multiple receivers did not affect the MUSIC peak am¬ 
plitudes, it does impact the effective SNR of the reconstructed signals, as evidenced by 
Fig. 22, where the two-element reconstruction is visibly noisier than the six-element one. 

Variations in SNR with multiple signals present causes additional distortions to the 
MUSIC DOA spectrum. The MUSIC spectra for two sinusoids of 1000 MHz and 
1005 MHz at 15° and 18° with SNRs of 10, 3, and —10 dB are provided in Fig. 23. As 
the SNR decreases, the ability to resolve the peaks disappears. Signals further separated in 
frequency have the same result, so the effect is likely solely due to the SNR and not to any 
signal interaction in the frequency domain. 

Next, the effect of signals with different SNRs was investigated. The MUSIC spectra for 
two sinusoids are presented in Fig. 24(a); the first is at 800 MHz and —5° and the second is 
at 1200 MHz and 15°. The SNR of the first is varied over 10, 0, and —10 dB, and the SNR 
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frequency (MHz) 

(a) Two-element array. 

Figure 22. Frequency spectrum of reconstructed 



frequency (MHz) 

(b) Six-element array. 

GHz sinusoid with received SNR of —10 dB. 


of the second is held constant at 10 dB. The first peak nearly disappears when that signal is 
20 dB less than the second. The large separation in DOA and frequency indicates that this 
is purely a SNR effect and is not caused by any frequency overlap. 


10 dB OdB -10 dB 



Figure 23. MUSIC spectrum effect of SNR on MUSIC spectrum DOA resolution with two 
sinusoids: 1000 MHz at 15° and 1005 MHz at 18°; average of Monte Carlo simulation with 100 
trials. 
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Angle off boresight (sinespace) Angle off boresight (sinespace) 


(a) Sinusoidal signals. (b) 10% bandwidth QPSK signals. 

Figure 24. MUSIC DOA spectra for signals with different SNRs. 

The same trial was re-run with QPSK signals with a bandwidth of 10% fc', the results 
are presented in Fig. 24(b). The peak amplitudes are reduced somewhat from the sinusoidal 
signals, as is expected based on results from Sections 3.3 and 3.4. In general, the drop in 
peak amplitude with increased SNR difference matches the results with sinusoidal signals. 
With wideband QPSK signals, there is an additional distortion of the MUSIC response of 
the high power signal. A possible reason for this distortion could be that since the array 
covariance matrix has very little information in the direction of the low-power signal, the 
low-power signal’s subspaces becomes indistinguishable from the noise subspace. Since 
the MUSIC algorithm is being run under the assumption that two signals are present, this 
distortion of the subspace results in an incorrect estimate of the noise eigenvectors and a 
degraded MUSIC response. 

3.6 Performance Analysis Review 

An analysis of the MUSIC algorithm has been performed with particular emphasis on the 
effects of signal frequency, bandwidth, DOA, and SNR. In most scenarios, MUSIC per¬ 
forms very well even with wideband signals. The use of a receive antenna array with sparse 
spacing and — 6 dB sidelobes likely contributes to some of the distortions and degradations 
in the MUSIC response, but the tradeoff in resolution over a uniformly spaced array may 
be worth the occasional inaccuracy, especially when averaged over multiple collections. 
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While this analysis was completed solely with the MUSIC algorithm, all of the sub¬ 
space algorithms are based on leveraging the same array covariance matrix eigen-structure 
and signal/noise subspace relationship. The same types of performance breakdowns as 
investigated above are expected with other basic narrowband subspace methods. 

In the next chapter, enhancements to the MUSIC algorithm are proposed to reduce the 
time required for computing a high-resolution result and to resolve ambiguous DOAs in 
excessively sparse arrays. 
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CHAPTER 4: 

MUSIC Performance Enhancements and Limitations 


4.1 Computational Performance 

One of the major limitations of the MUSIC algorithm is that the seareh over the angular 
visible range is eomputationally intensive. The angular step size utilized in that seareh 
contributes to DOA error, so a small angular step size is desirable. A small step size means 
that there will be many angular values for which to test the orthogonality via (2.21). 

4.1.1 Computational Complexity 

Given a collector with ap-element receiver array, the array covariance matrix will be p x p, 
and the eigenvectors will be of length p. The number of noise eigenvectors (E^) will de¬ 
pend on the number of signals; in general there will ht p — q noise eigenvectors. The 
steering vectors (a(0) will also be of length p. Neglecting any computations to derive the 
steering vector components and assuming that the steering vectors are already normalized, 
we only need to account for the matrix multiplications in the denominator of (2.21). Defin¬ 
ing complexity as the number of floating point operations and assuming that both addition 
and multiplication are complexity order one 0{\), the complexity of a multiplication of 
two matrices of sizes mx n and n x k is 0{mk{2n — 1)) ^ 0{2mnk). The approximate 
complexity of the MUSIC algorithm for a single DOA is then: 

0{MUSIC) =0 [a^(0)EwEw"a(0)] 

=0[{1 X p){p X {p - q)){{p - q) X p){p X 1)] 

=0 [2p{p -q) + 2p{p -q)+ 2{p - q)] 

=0{Ap^ - 4pq + 2{p-q)) 

^0{4p^ — 4pq) (4.1) 

The number of receive elements will generally be relatively small, especially for a sparse 
array designed for aerospace applications; therefore, the leading scalars and the second 
order terms are retained because they will still have a large impact on the result. As an 
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example, for the six-element array of (3.1) with a single signal present the eomplexity is 

~ 0 ( 120 ). 

This portion of the eomputational eomplexity is determined by the MUSIC algorithm; 
redueing it would require use of a different direetion finding algorithm. The number of 
diserete angles to seareh and the angular step size between them is at the diseretion of the 
system designer, however, and this provides an opportunity for improving eomputational 
performanee. In general, if an angular resolution of A is desired, the number of diserete 
angles over whieh to eonduet a one-dimensional MUSIC seareh in azimuth is m = a/A, 
where a is the angular width of the visible range. In the analyses of Chapter 3, a is 180°, 
or two in sinespaee. For a direet implementation of MUSIC then, the eomplexity order for 
determining the full DOA speetrum would be O (m (4p^ — Apq )). For praetieal appliea- 
tions, the number of diserete angles m will generally be mueh larger than the number of 
reeeivers and will be the primary driver of the number of eomputations required, espeeially 
when a two dimensional (azimuth and elevation) solution is desired. Therefore, redueing 
the number of diserete angles ean be of great advantage to improving proeessing speed. 

4.1.2 Reducing Search Complexity 

One method of redueing the number of diserete angles without also redueing the resolution 
would be to run the MUSIC algorithm more than onee. On the first run, a large step size 
eould be used whieh provides only a rough estimate of signal direetion. Those estimates 
ean then be used to re-run MUSIC with a mueh smaller step size over only the angles 
surrounding eaeh of the signal direetion estimates, with the result of obtaining a high reso¬ 
lution DOA estimate without the large m whieh would normally be required. The question 
arises then of how best to ehoose a eourse resolution step size. 

Sinee the ideal MUSIC DOA peak is infinite, a logieal deeision is to ehoose a step size 
whieh will provide a “good enough” peak height to distinguish signal DOA peaks from 
anomalous peaks whieh arise due to array sparsity in a linear array. Examination of the 
various MUSIC speetra plotted in Chapter 3, and espeeially that provided in Fig. 11(e) 
indieate that, at least for the array layout of (3.1), the anomalous peaks generally do not 
exeeed a 10 dB orthogonality. While this is far from an exhaustive evaluation and is only 
based on trials with three or fewer simultaneous signals, it is aeeeptable for the purposes of 
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this discussion. An appropriately separated signal peak may then be defined as having an 
orthogonality of at least 20 dB. A step size can be derived from this by measuring the width 
of a signal’s high resolution MUSIC response at a 20 dB orthogonality. The width and 
height of the MUSIC response of wideband signals can vary depending on the bandwidth 
and DOA but will never be narrower than the response for a sinusoidal signal because the 
minimum width is related to the array beamwidth, as indicated in the results Section 3.1.1. 
Therefore, the width of a sinusoidal signal at 20 dB orthogonality provides an acceptable 
metric for a coarse step size, Acrs- A fine resolution DOA can then be obtained by running 
(2.21) a second time with a fine step size over the range of [—Af„/2, Ac„/2). 

A portion of the MUSIC DOA spectrum in a fine and two coarse angular step size 
for two sinusoidal signals at —29° and —12° is provided in Fig. 25. The fine step size 
is 2 X 10^^ in sinespace (~ 0.01° at boresight), and the first coarse step size is 0.01 in 
sinespace (0.6° at boresight), which is the width of the fine step MUSIC response at 20 dB. 
The ratio of those step sizes is 50:1; while the exact advantage of using a two-step approach 
will depend on a and the number of signals present, a nearly 50 times speed improvement 
is highly desirable. As an example, when a = 2 in sinespace and two signals are present, 
the fine step MUSIC requires 10,000 executions of (2.21), and the two step approach would 
require only 300, which is 33 times faster with the same signal DOA resolution! A third 
step size of 0.025 is also provided in Fig. 25 which indicates that as the step size increases, 
a point is reached where the MUSIC peak locations may no longer be resolvable. 

It is also worth looking at the geolocation error induced by these step sizes. For an 
aircraft at 30,000 ft the horizon is ~180 nmi away. An angular step of 2 x 10^^ implies 
a location uncertainty of ±110 ft at that distance, while an angular step of 0.01 implies a 
location uncertainty of ±0.9 nmi. If many collections of the same signals from the same 
direction were able to be completed with varying step sizes or step size offsets, this error 
could be averaged out. The more realistic scenario is that only a few signal collections 
may be likely from different positions as the collector moves, a situation in which use of a 
coarse step size could result in a significant location bias error. 

While not a reduction in search complexity, it should be noted that the MUSIC algorithm 
is ideal for parallelization, which can reduce computation time when multiple computa¬ 
tional cores are available or when the process is being embedded in a field programmable 


47 




Angle from boresight (sinespace) 


Figure 25. MUSIC DOA spectrum with varying angular step size. 

gate array (FPGA). Since the MUSIC algorithm is a search over DOAs and the results for 
each discrete DOA are independent, each available core can be tasked with calculating a 
portion of the visible range. There is also no particular limit to the number of cores over 
which the process could be distributed aside from the number of discrete angles. 


4.2 Ambiguity Resolution 

In his paper on the MUSIC algorithm, Schmidt notes that there is an inability to resolve type 
I ambiguities, where a{0i) = a{ 62 ) and 0i ^ O 2 [2]. This is true for a stationary collector 
and emitter, but ambiguity resolution can be accomplished when the collector is moving 
at a known velocity and multiple collections of the same signal can be obtained. This is 
because the true DOA will always point in the direction of the emitter, but any ambiguities 
will tend to drift. 


To illustrate this ambiguity resolution ability, the array layout of (3.1) has been doubled 
so that the minimum spacing is now one wavelength: 


d 


0 1 3 7 12 20 


wavelengths. 


(4.2) 


This array spacing now has a grating lobe in addition to the main beam, which is illustrated 
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in the array radiation pattern plot of Fig. 26, where the main beam is at 10° and and the 
grating lobe is at 55.7°. This grating lobe manifests as a DOA ambiguity; there is no method 
of determining from a single collection which is the correct DOA. It should be noted that 
in sinespace the mainlobe/grating lobe separation will always be the same regardless of the 
mainlobe’s offset angle from broadside; this is not true when plotting the radiation pattern 
in degrees. For a minimum receiver spacing of one wavelength, the lobes will be separated 
by one in sinespace. 



Figure 26. Radiation pattern for array with main beam steered to 10° and grating lobe at 55.7°. 

To make an estimate of emitter location, all DOA vectors collected for a particular 
signal are utilized simultaneously. To determine DOA drift, however, only two need to 
be examined—the current and previous. By keeping track of the intersection points of 
DOA vectors from subsequent collections, it should easily be determinable that one set of 
crossing points remains relatively fixed while another set drifts. 

An illustration of this principle is provided in Fig. 27, with a depiction of 2-5 collections 
of two 0 dB SNR sinusoidal signals near fc. Each signal is a different color, and the two 
DOA vectors for each signal at each collection represent the two possible DOAs, only one 
of which is actually correct. The circles represent the DOA intersections of each subsequent 
collection. 
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(a) Two collections. 
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(c) Four collections. (d) Five collections. 

Figure 27. Ambiguity resolution by tracking of DOA intersection points. 


From examination of Fig. 27, it can be seen that the intersection point of one set of 
DOA vectors remain fairly constant, while the intersection points of the other set drift 
fairly significantly. The true DOAs are easily distinguishable by five collections, and a 
reasonable guess can be made with only three collections. 

The MUSIC DOA spectra for the previous example are provided in Fig. 28, and it can 
be seen that the spectrum over sinespace [—1,0) is an exact replica of that over sinespace 
[0,1). This is due to the one wavelength minimum receiver separation in the array and 
the sinespace one separation in antenna lobes. When an ambiguity is present in the array. 
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one no longer has to process the entire visible range since each signal will be repeated 
more than once in the MUSIC spectrum. This does not necessarily equate to halving the 
processing time; since the antenna beamwidth is now narrower, the MUSIC response is 
narrower, and a smaller step size may be desired to ensure an appropriate peak amplitude 
can still be determined. 



Figure 28. MUSIC DOA spectra for receiver array with single ambiguity. 

A metric for comparing the motion of intersections between DOA possibilities can be 
calculated by determining the vector between any two successive intersections, summing 
all of these vectors for each DOA, and determining the square of their length. This is also 
equivalent to just determining the magnitude of the vector from the first to the last collected 
intersection points. The metric for the true DOA should be much smaller than that for any 
ambiguous DOAs since the ambiguous DOA intersection points will travel further. The 
metric values for the prior example are provided in Table 7. By these metrics, for Signal 1 
(blue in Fig. 27) the negative DOA is correct, and for Signal 2 (red in Fig. 27) the positive 
DOA is correct, which matches the graphical interpretation. 

As shown in Chapter 3, high bandwidths and signals of different power can result in var¬ 
ious sorts of DOA estimation degradations, of which the severity is dependent also on the 
exact DOAs and frequencies present. These degradations lead to DOA estimation errors, 
which can lead to the intersection points of subsequent signals not to drift in an orderly fash- 
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Table 7. Ambiguity metric over multiple collections with sinusoidal signals. 




Collection 

Signal 

DOA 

3 

4 

5 

Signal 1 

DOA 1 

185 

391 

10 

DOA 2 

1408 

5869 

10646 

Signal 2 

DOA 1 

1093 

4244 

11191 

DOA 2 

152 

205 

6 


ion or to jitter around the true emitter position. A slightly more robust ambiguity metrie ean 
be ealeulated using the midpoint of two subsequent intersection points vice the intersection 
points themselves. One downside is that this requires at least four collections before mak¬ 
ing a judgment on which is the correct DOA. An illustration of the map view and MUSIC 
spectra for six collections with a 15 percent bandwidth LFM and QPSK signal, both at 
10 dB SNR, are provided in Fig. 29(a). Circles again indicate the intersection points and 
squares indicate the mid-position of subsequent intersection points. The ambiguity metrics 
calculated using both the DOA intersection points and the intersection point mid-positions 
are provided in Table 8, and a ratio of the metric for the true DOA divided by the metric for 
the ambiguous DOA is provided for easy comparison between the two methods. In general, 
the mid-position metric provides a better estimation and has less variability. 




(a) Map view. 


(b) MUSIC DOA spectra. 


Figure 29. Ambiguity resolution for receiver array with single ambiguity and wideband signals. 
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Table 8. Ambiguity metrics over multiple collections with wideband signals. 






Collection 


Signal 

Point type 

DOA 

3 

4 

5 

6 



DOA 1 

456 

6857 

9708 

8878 


Intersection 

DOA 2 

58 

47 

2334 

293 

Signal 1 


Ratio 

0.127 

0.007 

0.240 

0.033 


DOA 1 

- 

1714 

6402 

7349 


Mid-position 

DOA 2 

- 

12 

289 

146 



Ratio 

- 

0.007 

0.045 

0.020 



DOA 1 

I 

54 

624 

10 


Intersection 

DOA 2 

852 

2397 

1975 

123985 

Signal 2 


Ratio 

0.001 

0.023 

0.316 

0.000 


DOA I 

- 

14 

245 

119 


Mid-position 

DOA 2 

- 

599 

1032 

33683 



Ratio 

- 

0.023 

0.237 

0.004 


One issue that arises with arrays with ambiguities is that the true DOA of one signal 
could be co-incident with an ambiguous DOA of another signal. In this case, the MUSIC 
response will only have one peak in that location, which could make separating the two 
signals difficult or impossible. 


4.3 Limitations 

While some enhancements and understanding of the MUSIC algorithm are possible, there 
are some fairly significant limitations which make real-world implementation difficult be¬ 
yond the computational complexity discussed in Section 4.1.1. 

4.3.1 Steering Vector Estimation 

Equation (2.8) is an analytical estimate for a steering vector based on the frequency, an¬ 
gle, and receiver mounting locations, but for a real array the steering vectors will not be 
linear functions of frequency and angle because of mutual coupling of receive antenna el¬ 
ements and effects of the vehicle structure to which the elements are mounted. For a real 
implementation, the steering vectors will generally be determined via array calibration and 
stored as a lookup table [2]. It would be impractical to measure and store a lookup table 
with very small steps in both frequency and angle: for an array with a 180° visible range, 
a 1 GHz bandwidth, a frequency step of 100 Hz and an angular step of 0.1°, a lookup table 
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with 1.8 X 10^ entries would be required—and that is only for azimuth determination! 

Errors in reeeive element positioning and pointing are also unavoidable, and even disre¬ 
garding mutual eoupling, the radiation patterns on multiple elements will not be identieal. 
All of these faetors eombine to ereate errors in one’s knowledge of the steering veetors 
whieh will result in additional DOA estimation errors. 

4.3.2 Signal Number Estimation 

One limitation that has been ignored in this paper is in estimating the number of signals 
q present during any eolleetion period. Knowing this is essential to determining whieh 
of the array covarianee matrix eigenveetors are signal eigenveetors and whieh are noise 
eigenveetors, beeause only the noise eigenveetors are used in the MUSIC algorithm of 
(2.21). This limitation is not speeifie to the MUSIC algorithm but exists for all subspaee 
algorithms. 

In Seetion 2.4, it is noted that the p — q smallest eigenvalues should all equal but 
in praetiee eigenvalues will fall along a continuum and it can be difficult to distinguish the 
smallest signal eigenvalue from the largest noise eigenvalue. There are numerous items 
in the literature which develop and analyze statistical likelihood ratio tests which attempt 
to solve this problem. Wax and Kailath wrote one of the seminal papers on the subject 
with a development and analysis of estimators based on two information theoretic criteria 
termed Akaike information criterion (AIC) and minimum descriptive length (MDL) [20]. 
Hill and Pickholtz later provided a Monte Carlo-based calculation approach to compare the 
effectiveness of a number of different estimator types [21], and more recently Nadakuditi 
and Edelman developed a new method simply referred to as the “new” method [22]. 

These methods are all functions of the number of receivers, the number of samples, and 
the eigenvalues of the array covariance matrix. Code to implement two of these methods. 
Wax’s MDL and Nadakuditi and Edelman’s “new” method, is contained in Appendix B. Eor 
the analyses in this paper, however, both estimators nearly always returned the maximum 
possible number of signals (five in the case of a six-element array) and were almost never 
correct. Due to time constraints further research on the reasons for this were not attempted. 
Eor this reason the estimated number of signals was set to equal the known number of 
signals for all implementation of the MUSIC algorithm employed in the examples and 
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analyses of this paper. 


4.3.3 Signal Matching 

Localization of a particular signal over multiple collections requires that the signals col¬ 
lected at different times can be identified as being the same signal. With MUSIC and other 
subspace algorithms, if more than one signal is present at one time the incident signals can 
only be estimated by reconstruction (see Section 2.4.2) using the estimated DOAs. Even 
with sinusoidal signals, reconstruction is only conducted using a limited number of sam¬ 
ples, and two collections will not match exactly. With radar pulses it may be particularly 
difficult to tell the difference between two emitters operating at the same frequency with 
such a limited sample set. Additionally, the blending that can occur with wideband signals 
(see Section 3.3) creates a severe impediment to accurate signal matching since the effect 
will be different for each combination of signals and DOAs. 

4.4 Performance Enhancements and Limitations Review 

In this chapter, a method of reducing the number of computations required to obtain a 
high-resolution DOA estimate was proposed; this method utilizes a two-step execution of 
the MUSIC algorithm to first obtain a coarse DOA spectrum over the entire visible range 
and then to calculate a fine DOA spectrum over very limited portions in the vicinity of the 
signal DOAs. A method of resolving ambiguities for moving receiver arrays with minimum 
receiver spacing > A/2 was also proposed, and two metrics were derived which can provide 
an estimate of the true DOA within three or four collections. Finally, some of the limitations 
inherent in the MUSIC algorithm and in this research were delineated. 
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CHAPTER 5: 
Conclusion 


5.1 Results 

This thesis research was conducted with three main goals: to implement the MUSIC sub¬ 
space direction finding algorithm and geolocation in MATLAB; to conduct a high-level 
performance analysis of the algorithm with particular emphasis placed on the effects of 
frequency, bandwidth, DOA, and SNR and to devise some enhancements to the process, if 
possible. 

5.1.1 Performance Analysis 

In the context of this research, performance is defined as the ability to resolve a signal’s 
DOA and the accuracy of that DOA. It has been shown that signals not at the array center 
frequency will have a shifted response in the MUSIC DOA spectrum, but the frequency 
alone does not result in a degradation of the response. It has also been shown that the 
weighted mean frequency, described in Section 3.2, provides an acceptable method of es¬ 
timating and correcting the MUSIC DOA shift. The accuracy of the subsequent DOA 
determination is then based on the accuracy of the reconstructed signal’s frequency. For 
wideband signals, it was shown in Section 3.3 that signal blending can occur during re¬ 
construction; this blending will corrupt the frequency estimation and the DOA estimation. 
Over multiple collections there may be situations where a signal-of-interest is blended with 
different signals on each collection; it is expected that this will result in a colored noise 
effect on the signal DOA estimates and the subsequent geolocation. 

Besides the effect on signal reconstruction, it was demonstrated that the MUSIC spec¬ 
trum peak amplitudes decrease roughly linearly with increased signal bandwidth. Because 
a peak still exists even at extremely high bandwidths of 40% fc, this reduction in ampli¬ 
tude might at first seem to be immaterial; however, when combined with an offset in DOA 
from array boresight, the amplitude reductions can be significant. Additionally, if multi¬ 
ple signals are present there will be a reduced ability to resolve two closely spaced signals 
and there will be a higher likelihood of distortion due to the presence of multiple signals. 
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An example of both of these effects was presented at the end of Section 3.3. While there 
are an infinite number of possible signal combinations and an exhaustive analysis was not 
conducted, this example indicates that 5-10 percent signal bandwidth is the effective limit 
for the baseline MUSIC algorithm to maintain signal resolution and avoid DOA distortion. 
The evaluation of the effect of DOA on the MUSIC spectrum in Section 3.4 also indicates 
that even with a bandwidth of < 10% fc, MUSIC is likely to only be truly effective out to 
~ 30° off boresight. This result is interesting also in that an array is likely to be comprised 
of receiver elements which have some directivity and do not have a perfectly semicircular 
pattern in azimuth. It may make sense to choose receive elements with beamwidths which 
effectively match the DOA degradation for the expected bandwidth of received signals. 
For the situation described above, receivers with an approximate gain of 9 dB and a 3 dB 
beamwidth of ~ 60° could be utilized. If one is willing to sacrifice the ability to detect nar¬ 
rowband signals beyond 30° off boresight, a computational benefit can also be realized in 
reducing the MUSIC search to only ±30°. Since 30° is 0.5 in sinespace, this is effectively 
a halving of the number of computations required as compared to a search over ±90°. 

An analysis of the effects of signal SNR was presented in Section 3.5, and it was shown 
that MUSIC is still fairly accurate even down to —10 dB SNR, though the peak amplitude 
is reduced. Monte Carlo simulations indicate that the effect of noise is unbiased over 
multiple collections, but that the accuracy of a particular collection DOA is degraded to 
a standard deviation of approximately a quarter of a degree at —10 dB SNR. Additional 
effects are discernible when multiple signals are present. First, as SNR decreases, there is a 
reduction in the angular resolution between source signals. Second, an investigation of the 
MUSIC response with overlapping signals of different power indicates that, especially with 
wideband signals, there is an effective limit of 10 dB difference between signal powers in 
order to be able to resolve the presence of the weaker signal in the MUSIC DOA spectrum. 

These results will prove useful to system designers considering implementation of the 
MUSIC algorithm. If the system will be used in an environment where signal overlap is ex¬ 
pected, source signals have a 5% bandwidth or less, and signal powers do not vary greatly, 
then the basic MUSIC algorithm should perform quite well. Between 5-10% bandwidth 
and depending on the number of concurrently overlapping signals, the basic MUSIC al¬ 
gorithm may not provide the desired accuracy, and as bandwidths of 20% or greater are 
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approached an alternate method which can correct for the bandwidth effects is likely re¬ 
quired. 

5.1.2 Enhancements 

The first performance enhancement developed in this thesis is a method of reducing com¬ 
putational complexity through the use of a two-stage DOA determination. In the first stage 
a coarse angular step size is utilized for the MUSIC search; the step size should be cho¬ 
sen to ensure a minimum peak amplitude for signal detection. The second step runs the 
exact same MUSIC calculation, but this time with a fine step size which should be chosen 
based on desired geolocation precision. The second stage is only run in the vicinity of the 
signal peaks which were extracted from the first stage. In this manner a high resolution 
can be obtained with significantly fewer calculations than if the fine step size were to be 
utilized over the entire visible range. The example developed in Section 4.1 resulted in 
an approximately 30 times speed improvement. This improvement enables the real-time 
implementation of the MUSIC algorithm in less powerful computing hardware than might 
otherwise be required or available. 

A geolocation enhancement was also introduced in Section 4.2. Generally, arrays with 
ambiguities, or grating lobes, are undesirable for direction finding applications. The mo¬ 
tion of an airborne or spaceborne array provides a method of resolving those ambiguities, 
however, under the assumption that the same signal can be received and correlated over 
multiple collection durations. This is possible because the true DOA will always point to¬ 
wards the transmitter, but the false DOA (s) will drift and not all cross in the vicinity of a 
single point. A DOA drift metric is also derived which enables automatic estimation of the 
correct DOA with a minimum number of collections. This enhancement enables the con¬ 
struction of arrays where receiver size or other constraints to the minimum receiver spacing 
result in the presence of unavoidable grating lobes. 

5.2 Future Work 

There are many areas in which this thesis work could be expanded. The most immediate 
research need would be in an examination of some of the more recent and more advanced 
subspace methods which could reduce computational complexity and reduce the impact 
of signal bandwidth on the determination of signal DOA. Friedlander and Weiss’s work in 
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using virtual arrays [23] should enable a reduction in computational complexity by enabling 
the use of root-MUSIC [7], but the accuracy of estimating interpolated signal phases for a 
sparse array would need to be closely examined, as would the complexities arising from 
the use of multiple virtual arrays over different portions of the visible range. There have 
been multiple papers written since Schmidt’s seminal paper on MUSIC [2] concerning 
methods to improve the subspace DOA resolution and minimize interaction effects for 
wideband signals; a primary candidate for future research would be the method developed 
by Friedlander and Weiss which builds on top of their prior work on virtual arrays [11]. 
A comparison to the methods of Doran, Doran, and Weiss [10] as well as other promising 
methods would also prove fruitful. 

There are also multiple possible avenues of directly expanding on the work done in 
this thesis. One area would be in an evaluation of the impact of collection duration on 
the MUSIC spectrum; at a minimum an improvement in DOA resolution at low SNRs is 
expected, but there may be other advantages and disadvantages as well. The implementa¬ 
tion of more realistic signals would confirm that methodology employed by this thesis of 
using representative test signals was appropriate; this should include more realistic signal 
durations and power levels. Examples could include phase-coded radar signals such as 
maximal-length and Frank polyphase codes or communications signals such as 16-QAM 
using Nyquist pulse shapes. Finally, the impact on DOA estimation and signal reconstruc¬ 
tion should be evaluated with irregularly timed signals which break across collections. 

Secondly, a deeper investigation into the methods of signal number estimation and their 
limitations could resolve the reasons why this thesis research was unable to obtain correct 
estimations. Being able to make an estimation with at least a reasonable expectation of 
accuracy is an essential piece of realizing an operational employable subspace-based di¬ 
rection finding system. This research should begin with the papers in [21] and [22] and 
include any other relevant work in the area. At a minimum, the impact of signal type, 
bandwidth, SNR, and array layout on the ability to estimate the number of signal should be 
investigated. 

A third research area would be in methods of improving the signal reconstruction to 
avoid the blending which can occur with wideband signals, as demonstrated in Section 3.3. 
An analysis of the ability to correlate separate collections as the same source signal in the 
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presence of blending could also prove useful. This is of particular concern in high-density 
signal environments where it can be especially difficult to separate and geolocation signals. 

A final possible field of study would be in applying the model error analysis of Swindle- 
hurst and Kailath [13] for non-linear arrays. Additionally, their weighting scheme for 
MUSIC could be evaluated to see how it performs for the same radar-type signals utilized 
in this thesis over a range of frequencies and bandwidths. 

5.3 Conclusion 

It is hoped that the results of this thesis, which are a first step towards deeper investigations 
into the subspace direction finding algorithms, will enable an expansion of the geolocation 
research which has been on-going at the Naval Postgraduate School to include more anal¬ 
ysis beyond the more traditional two-receive element geolocation methods. The ability to 
accurately utilize MUSIC in moderate to low bandwidth environments, despite it being de¬ 
rived specifically for sinusoidal signals, was also demonstrated. Additionally, the density 
of RF emitters is increasing around the world, and with that comes an increase in the like¬ 
lihood of receiving overlapping signals. Methods of accurately resolving and geolocation 
those overlapping signals are likely to be of interest to the military in the future, which 
means that now is the time for research into subspace direction finding and other promising 
methods of multiple signal geolocation. 
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APPENDIX A: 

Subspace Processing Specifics 


A.l Time-Domain Steering Matrix with Wideband 
Signals 


The time-domain solution derived in Section 2.4, and specifically the definition of the steer¬ 
ing vectors as (2.8), is only an exact relationship for sinusoidal signals. To show the true 
nature of the steering matrix, it is useful to derive an analytical solution for the representa¬ 
tive wideband signals utilized in this thesis. 

Let (t) be a LFM pulse at center frequency fd, bandwidth B, and duration T, and let 
52(t) be a phase-coded signal of frequency fd and phase sequence 

(A.l) 

(A.2) 


Additionally, assume that the two signals arrive at time zero at the origin of a three-element 
array. There will be no phase shift at the origin, but the phase shift of signal st at a receive 
element offset dj from the origin will be in accordance with (2.1). For this derivation, it 
makes more sense to express the phase in terms of frequency, where c is the speed of light: 


(j)ij = iTtfi—sinei. 


(A.3) 


The DOA 6i and the receiver offset from the origin dj are both straightforward, but the 
frequency f is not. The instantaneous frequency of a signal is the derivative of its phase in 
cycles; for the two signals in (A.l) and (A.2) the phase in cycles is 


h {t) = (^fci - f) ^ (A.4) 

hit)=fc2t+^Wit)- (A.5) 
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The instantaneous frequency of each signal is then 


flit) = + 

(A.6) 


(A.7) 


It is clear that the instantaneous frequency of each signal is time-dependent, potentially 
non-linear, and will change depending on the specific frequency/bandwidth/phase parame¬ 
ters of each signal. Substituting (A.6) and (A.7) into (A.3) will result in a time-dependent 
steering vector. The constant steering vector of (2.8) is then simply an approximation 
which may only truly be accurate over a very short portion of the collection. It is from this 
inaccuracy that most of the MUSIC degradations explored in Chapter 3 arise. 


A.2 Subspace Methods in the Frequency Domain 

The array matrix relationship of (2.6) can also be developed in the frequency domain; doing 
this offers some computational advantages over the time domain method. At the origin, the 
received signal is the sum of all q incident signals, just as in the time domain, with no phase 
shifts: 

Xoif) = Y,Siif). (A.8) 

i=i 

where Si{f) is the Fourier transform of the zth signal. For any particular source signal 
Si{f), there will be a phase shift in the received signal at a receiver offset from the origin, 
and that phase shift will be the same as it is in the time domain. In Section A.l, it was 
shown that the frequency component of the phase shift was time dependent for wideband 
signals and had to be derived separately for each signal. In the frequency domain, however, 
assuming that the DFT was utilized to obtain the frequency spectrum, the frequency term of 
the phase shift will simply be the vector of frequency bins from the DFT. This simpler and 
always linear expression of the steering vectors for wideband signals is the first advantage 
of conducting the subspace algorithms in the frequency domain. A second advantage is 
that even though there are the same number of frequency bins as time domain samples, it 
is relatively easy in the frequency domain to reduce the number of samples for which the 
algorithm needs to be run. The most obvious way to reduce the number of samples in the 
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calculation of the covariance matrix in the frequency domain is by using only the positive 
frequency bins. Since the subspaee algorithms are run on the analytic signal, as long as the 
algorithm is executed at RF or IF the negative frequency portion of the real signal frequeney 
spectrum can be disearded sinee it is simply a duplieate of the positive frequeney portion. 
The number of eovariance samples can be further redueed if one is only eoncemed about 
signals over a limited frequeney band. The main potential issue with reducing the number 
of reeeived signal samples in the frequeney domain is that there are then fewer samples 
over whieh to average out the noise, and the noise in any particular frequency bin can have 
a larger impaet on degrading the signal and noise subspaee relationship. 

In fact, since the Fourier transform is eomplex and can be expressed in magnitude and 
phase eomponents, only the phase portion of the frequeney speetrum is neeessary to eom- 
pute the array eovariance matrix and determine signal DOAs. Sinee the reeeived signal 
at an element offset from the origin is simply an additional phase shift, the magnitude of 
the reeeived frequeney speetrum is exaetly the same at all receivers, negleeting noise and 
very small differenees in distanee from the souree and assuming that all receivers have an 
identieal radiation pattern. 

Working in the frequeney domain also means that if any signal frequency corrections 
are required (see Section 3.2), all of the required information to estimate signal center 
frequeneies is already present and one does not have to switeh baek and forth between the 
time and frequeney domains. 

A.3 Subspace Method Processing at an Intermediate 
Frequency 

Subspaee method direetion finding proeessing ean be eompleted at an IF with the same 
results as if the processing were eompleted at RF, disregarding any non-linearities resulting 
from the mixing proeess and assuming that the reeeived signal at all reeeiver elements is 
mixed down to the exact same IF prior to sampling. For the purposes of this diseussion, fr 
is the array center frequeney at RF, /„ is the mixer frequeney, and fi = fr — fm is the center 
frequeney of the IF reeeived signal. 

When using the MUSIC algorithm at IF, while the array eovariance matrix (2.12) is 
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computed using the IF signals, the steering vectors used in the DOA search of (2.21) must 
be computed using the RF array center frequency fr. This same steering vector construction 
with fr is used during signal reconstruction and frequency estimation. The resulting signal 
estimated frequency in this case will be at IF, but any frequency correction (see Section 
3.2) must be completed as it would be for the RF signal by including an adjustment for the 
mixer frequency. Therefore (3.5) must be re-written as 


digest = arcsin 


-—sin0;',Mt/5/c ) • 

Ji,est “T Jm / 


(A.9) 
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APPENDIX B: 
MATLAB Code 


A working example of the MUSIC direetion finding algorithm, along with test signal gen¬ 
eration, estimated signal reeonstruetion, geoloeation, and plotting of results is provided 
below. The enhaneements diseussed in Chapter 4 are also ineluded in the primary seript. 
A supplemental eleetronie eopy of this eode ean be obtained by eontaeting the Naval Post¬ 
graduate Sehool’s Dudley Knox Library. The eleetronie version of the eode also ineludes 
a seript whieh generates the example and visualization eontained in Chapter 2 and a seript 
and assoeiated signal definition files whieh generate the results of Chapter 3. 

B.l MUSIC Direction Finding and Geolocation Script 


%% MUSIC direction finding algorithm and geolocation 

o, 

o 

% Purpose: 

% This script generates test signals which are defined by a specified 
% file and runs the MUSIC algorithm in the frequency domain in azimuth 
% only based on parameters sepcified in the initilization section below. 
% Signal reconstruction based on the derived DOAs is also conducted, and 
% if more than three collections are commanded an estimate is made of 
% the emitter location. Data is presented in plots for easier analysis. 

o 

0 

% Custom functions called: 

% sig_config_xxxx - not a function, but a script file which contains 
% specifications for one or more signals as well as the array center 
% frequency, sample frequency, and signal duration 
% interferometry_sig_gen - function which generates phase shifted 
% signals at specified receiver locations 

% subspace_calc - function which generates the array covariance matrix 
% and the eigendecomposition of said matrix. Also estimates number of 
% signals if desired 

% music_calc - computes MUSIC spectrum 

% geolocate2 - 2D least squares estimate of emitter locations based on 
% collected DOA vectors 
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% pt_motion - metric for ambiguity resolution; measures the motion of a 
% test point (currently intersection mid-positions) over multiple 
% collections 

o, 

o 

% History: 

% DATE (YYMMDD) AUTHOR/EDITOR EMAIL NOTE 

% 141203 Chris Straessle gcstraes@nps.edu Original thesis code 

o, 

o 

% NOTE 

% This is a full implementation of the code used for this thesis, but 
% it does not directly generate the figures contained in Chap 3 of the 
% thesis. Those were obtained through rearrangements and subplotting in 
% order to highlight specific results; the script used for that is 
% attached electronically to the NPS library's electronic version of 
% this thesis. The code used to generate the example and visualization 
% in Chap 2 of the thesis is also attached electronically. 

o, 

o 

9-9'9-9-9-9-9-9'9-9-9'2-9'9-9-9-9-9-9'9'Q-9'9-9'9-9-2-9'9-9'9'9-9'9-9'2-9-9-9'9-9'9'Q-9'9-9'9-9-9-9'9-Q-9'9-9'9-9'9-9-9-9'9-9-9'Q-9'9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%% Initalization - User options 


clear all 
close all 
%clc 

% options are 'uniform' or 'sparse' 

% 0 = no ambiguities, integer = # of ambiguities 
% Turn noise on or off. 1 = on, 0 = off 
% 1 = save figures, 0 = do not save figures 


arraytype = 'sparse'; 
num_ambig = 0; 
noiseon = 0; 
save_figs = 0; 


% compute signal number estimation. Options are 'none', 'NEW, or 'MDL' 
sig_num_est = 'none'; 


sig_config_example; 
selected_sigs = [1 2]; 


% load signal configurations from file 
% choose signals to model 


% range of frequencies for 
fmin = recParam.centerfreq 
fmax = recParam.centerfreq 


array (Hz) 
- 500e6; 

+ 500e6; 


68 




62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 
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nuin_collects = 5; % number of collection runs 

%% Additional Setup 
global c 

c = 3e8; % speed of light (m/s) 

% Angles at which MUSIC spectrum will be computed {this is main driver 
% of processing time). Two-stage method is being utilized. 
sine_step_crs = 0.001; % Coarse step size (sinespace) 

sinespace = -1:sine_step_crs:l-sine_step_crs; % discrete angles 
angles = asind(sinespace); % discrete angles (deg) 

sine_step_fine = 2e-4; % fine step size (sinespace) 

% initializations 
shift_req = 0; 

% custom color order definition 


[ 0 

0 

1.0000; 

1.0000 

0 

0; 

0 

0.5000 

0; 

0.7500 

0 

0.7500; 

0 

0 

0; 

0.7500 

0.7500 

0; 

0.2500 

0.2500 

0.2500] 


%% Array Configuration 

% Receive element locations; array spacing in the specified vector is 
% center frequency wavelengths 
switch arraytype 

case 'uniform' % half wavelength spacing 
rec_el_NL_wav = ([0 1 2 3 4 5]'/2); 
case 'sparse' 

% sparse spacing (minimum half wavelength spacing, no 
% ambiguities, consistent sidelobe height of ~ -6 dB). Pattern 
% is to place additional emitters at distances which equate to 
% an unbroken string of integer half wavelengths when all 


in 
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% elements are considered relative to each other. For the 
% current 6-element array, the half wavelength separations 
% covered are: 1-9, 11-13, 17, 19, 20. The first skipped integer 
% is 10, so the next element should be added 10 half wavelengths 
% after the last one. Adding a 7th element at d=30/2 would add 
% the separations: 10, 18, 23, 29, 30. Array can continue to be 
% expanded in this manner. 
rec_el_NL_wav = ([0 1 3 7 12 20]'/2); 
otherwise 

error (’ Incorrect array type specified') 

end 

% adjust array spacing for desired number of ambiguities 
num_lobes = num_ambig + 1; 

rec_el_NL_wav = rec_el_NL_wav * num_lobes; 

% Receive element locations in meters 

rec_el_NL = rec_el_NL_wav9rc/recParam. centerf req; 

p = length{rec_el_NL); % number of receive elements 

% Steering vector definition. Defined as a function for later 
% flexibility 

Af = @{f,ss) exp (-1 j*29cpi*f9rrec_el_NL9rss/c) ; 

% Determine number of ambiguitities at center frequency 

phi = 2*pi*rec_el_NL_wav*sinespace; % matrix of receiver phases by angle 
V = exp{1j*phi)/sqrt(p); % complex phases, corrected for array size 

% normalized complex antenna pattern at specified angle, shift to 
% 360*3/7 deg (increases peak finding accuracy); using a large enough 
% prime number denominator (7) minimizes the chances of a peak wrapping 
% around at H—90 deg (which is difficult to count properly); the 
% numerator (3) servers to more or less recenter the main beam, 
pattern = v(:,ceil(length(sinespace)*3/7))'*v; 
pattern_mag = abs(pattern); 

% ambiguity is being defined as > 90% of main beam magnitude. The MUSIC 
% response for grating lobes below 90% is small. This step truncates 
% the data at 0.9 
pattern_mag(pattern_mag<=.9)=0; 
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[pks. Iocs] = findpeaks(pattern_mag); 

nuin_lobes = length (pks); % count number of peaks remaining 

%% Collector details 

coll_pos = [0; 0]; % initial collector position 

coll_vel = [300; 0]; % nmi/hr 

coll_vel = coll_vel/3600; % nmi/s 

collect_interval = 120; % collection interval, in seconds 

% rotation matrix to normalize relative DOAs 
% (not currently utilized) 

R = 1/norm{coll_vel)*[coll_vel{1),-coll_vel(2);coll_vel(2),coll_vel(1)]; 
%% Collection Loop 
for n = 1:num_collects 

% collection time (s); 1st collection defined as t = 0 
t_col(n) = (n-1)*collect_interval; 

% collector position 

coll_pos ( : , n) = coll_pos { : , 1) + coll_vel9ct_col (n) ; 

% signal generation 

[sigs, X, e_pos, sig_names, doas(:,n), N, t] =... 

interferometry_sig_gen(sigParam, recParam, selected_sigs, ... 

coll_pos(:,n), coll_vel, rec_el_NL, noiseon) ; 

q = length(selected_sigs); % number of source signals 

Xf = fft(X,[],2); % convert to freq domain 

df = ((0:(N-1)).*(recParam.sampfreq/N))'; % FFT frequency vector 

% Subspace determination 

[q_est, EigVec_s, EigVec_n, EigVal] = subspace_calc(Xf, ... 
sig_num_est, q, n) ; 

%% MUSIC 
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% Run MUSIC with coarse steering vectors 
music_spectrum{n,:) =... 

music_caic{Af{recParam.centerfreq,sinespace), EigVec_n); 

% define the number of desired peaks to equal the number of 
% estimated signals times the number of ambiguities 

% NOTE - This assume that each true DOA has a grating iobe, which is 
% only exactly true for integer muitipies of minimum wavelength 
% spacing 

num_peaks = q_est*num_iobes; 

% extract peaks from MUSIC plot in descending order 

% NOTE - high bandwidth signals can lead to an anomoious peak being 
% chosen over a true peak. 

[pks, iocs] = findpeaks{real{music_spectrum(n, sortdescend ’) ; 

% re-sort q_est tallest peaks into DOA order {negative to positive ), 
% store as coarse angle vector 

ambig_iocs_crs_ss{n} = sinespace(sort(iocs(I:num_peaks))); %sinespace 
ambig_iocs_crs_ang{n} = angles(sort{iocs(I:num_peaks))); % degrees 

% re-run music calculation with higher precision steering vectors 
for m = i:num_peaks; 

% angle around which to conduct higher resolution search 
peak_ang = ambig_iocs_crs_ss{n}(m); 

% define fine resolution segment of visible range 
sine_sub_ang = (peak_ang-sine_step_crs/2): ... 

sine_step_fine:{peak_ang+sine_step_crs/2); 

% re-run MUSIC 
music_hires(n,:) =... 

music_caic(Af(recParam.centerfreq,sine_sub_ang),EigVec_n); 

% determine appropriate higher resolution DOAs and store as 
% precise angle vector 

% NOTE - currently written to only Took for one peak 
[-, indx] = max(real{music_hires(n, :))) ; 

ambig_iocs_p_ss{n}(m) = peak_ang - sine_step_crs/2 +... 
indx*sine_step_fine; 

end 

% convert to degrees 
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ambig_locs_p_ang{n} = asind (ainbig_locs_p_ss {n}) ; 

% fit DOA estimate into an ambiguities x signals matrix 
doa_est_ss{n} = reshape{ambig_locs_p_ss{n},q_est,num_lobes; 

%% frequency Domain Signal Reconstruction 

% use array center frequency and derived DOAs to solve for an 
% estimate of the source signals 

SIG_est{n} = pinv(Af(recParam.centerfreq,doa_est_ss{n}(1,:)))*Xf; 

% convert to time domain 

sig_est{n} = ifft(SIG_est{n}, [] , 2) ; 

%% Frequency Estimation 

% create normalized PSD for each estimated source signal 
SIGPWR_est_norm = (diag 

max(abs(SIG_est{n}(:,l:ceil(N/2))),[],2))+.. . 
abs(SIG_est{n}(:,l:ceil{N/2)))).^2; 

% truncate lowest 10% of response {minimizes noise and effects of 
% signals not being centered in search bandwidth) 
SlGPWR_est_norm(SIGPWR_est_norm<.1) =0; 

% estimate signal center frequency 

f_sig_est{n} = ((SIGPWR_est_norm*df(l:ceil(N/2) ))./... 
sum(SIGPWR_est_norm,2)).'; 

%% MUSIC adjustments 

% fix DOAs based on incorrect frequency assumption, uses fft 
% estimated center freq from above 
doa_est_ss_ff(:,:,n) = recParam.centerfreq ./.. . 

repmat(f_sig_est{n},num_lobes,1).*doa_est_ss{n}; 

% Determine direction vector to emitter 

% TBD: and correct for movement not in the +x direction 
% w = {# coll}(sig #, ambig #, (sin(DOA)|cos(DOA)) 

% Note that (sin(DOA),cos(DOA)) is an (x,y) unit vector in the 
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% direction of the emitter (assuming +x motion) 
w{n} = doa_est_ss_ff(:, : , n) ; % x components 

w{n}{:,:,2) = sqrt(1-w{n1).^2 ) ; % y components 

%% geolocation determination 

% NOTE - difficult to functionalize as currently written) 

if n > 1 

% intersection point of two DOAs 
intersect_pos{n} = intersect_pt{coll_pos,w); 
if n > 2 

% if no ambiguities, solve geolocation problem 
if num_lobes == 1; 

doa_est_ss_ff_na = reshape{doa_est_ss_ff,q_est,n) ; 
e_pos_est = geolocate2{coll_pos, doa_est_ss_ff_na); 

else 

% if ambiguities, attempt to determine correct DOA 

mid_pos{n} = {intersect_pos{n} + intersect_pos{n-1})./2; 
if n > 3 

% note: length(sum(distance vectors) has more dynamic 
% range than sum(length(distance vectors). Also, use of 
% mid-position has better results than individual 
% crossing points (though the latter can be calculated 
% one iteration sooner) 

sq_mid_motion{n} = pt_motion(mid_pos); 

% min values are likely correct DOA (not ambiguity) 
[~, indx] = min(sq_mid_motion{n}, [ ],1); 
for m = 1:q_est 
if n == 4 

% go back and remove ambiguous DOAs from 
% first 3 collections 
for k = 1:3 

doa_est_ss_ff_af(k, m) =... 
w{k}(indx(m),m,1); 

end 

end 

% remove ambiguous DOA from current collection 
doa_est_ss_ff_af(n,m) = w{n}(indx(m),m,1); 

end 
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end 

end 

end 

end 

end 

%% plots 

% determine if signals 'wrap' around ^—90 deg; shift DOAs if necessary 
% (note - cumulative if there's more than one wrap) 
for n = 2:num_collects; 

if ambig_locs_crs_ss{n}(1) > ambig_locs_crs_ss{n-1}(1) 

shift_req = shift_req + 1; 

end 

end 

if shift_req ~= 0; 

ambig_locs_crs_ss{n} = circshift(ambig_locs_crs_ss{n}, [0,shift_req]); 

end 

% MUSIC plot 
fignum{l) = figure{); 

plot (sinespace,10*logl0(abs{music_spectrum))) 
grid on 

xlabel {' Angle from boresight (sinespace) ', 'fontsize' ,14) 
ylabel (' Orthogonality (dB) ', 'fontsize ' ,14) 
set(gca, 'xtick' ,-1:.2:1) 
if num_collects > 1 

legend( 'c_l' , 'c_2' , 'c_3' , 'c_4' , 'c_5' , 'c_6' , 'c_7' , 'c_8' ) 

end 

title('MUSIC Specta', 'fontsize ', 14) 

% Signal FFT plot 
fignum(5) = figure(); 

SIG = fft(sigs.'); 
stem(df/le6,abs(SIG)) 

xlim([fmin fmax]/le6);set(gca, 'xtick' , (fmin:le8:fmax)/le6) ; 


Estimate emitter position using estimated correct 
DOA 

.pos_est = geolocate2(coll_pos, doa_est_ss_ff_af); 
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grid on 

xlabel{ 'frequency (MHz)’font size', 14) 
ylabel{ 'Volts (V)' , 'fontsize' ,14) 
legend {sig_naines) 

title (' Transmitted signals', 'fontsize' ,14) 

% plot locations and estimated DOA lines 

% signals are different colors and ambiguities are the same color 
cnum = num_collects; 
fignum{2) = figure{); 

plot(e_pos(1,:),e_pos(2,:), 'k+' , 'markersize' ,10) 
hold on 

plot(coll_pos(1,1:cnum),coll_pos(2,1:cnum), 'ko' , 'markersize' ,8) 

xlim([-250 250]) 

ylim([-25 300] ) 

xl(l,l,:) = coll_pos(1, :) ; 

xl = repmat(xl, [1 q_est 1]); 

= coll_pos(2,:); 
yl = repmat(yl, [1 q_est 1]); 
for n = 1:cnum 

for m = l:num_lobes; 

% wO = {# coll}(# ambig, # sig, sin(DOA)|cos(DOA)) 
wO{n} = doa_est_ss{n}; 

w0{n}(:,:,2) = sqrt (1-wO {n }(:,:, 1).''2 ) ; 

endpt{n}(m,1:size(wO{n},2),:) = wO{n}(m,:,:)*300 + ... 

repmat(reshape(coll_pos(:,n),l,l,2), [1,q_est,1]); 

% plot initial MUSIC DOAs 

line([repmat(coll_pos(1,n),1,q_est); endpt{n}(m,:,!)],... 

[repmat(coll_pos(2,n),1,q_est);endpt{n}(m,:,2)],... 

'linestyle' , '-' ) 
if n > 1 && num_lobes ~=1; 

plot(intersect_pos{n}(m,:,1),intersect_pos{n}(m,:,2),... 

'o','color', colors(n-l,:)); 
if n > 2 

plot(mid_pos{n}(m,:,1),mid_pos{n}(m,:,2), 'ks' ); 

end 

end 

end 

end 
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xlabel{ 'In-track Range (nmi )' , 'fontsize' ,14) 
ylabel {' Cross-track Range (nmi)'fontsize' ,14) 
title( 'Map View’ , 'fontsize * ,14) 

% plot locations and estimated DOA lines for frequency fixed DOAs 
% signals are different colors and ambiguities are the same color 
fignum{3) = figure{); 

plot(e_pos(1,:),e_pos(2,:), 'k+' , 'markers!ze' ,10) 
hold on 

plot(coll_pos(1, :),coll_pos{2, :), 'ko’ , 'markersize' , 8) 
xlim([-250 250]) 
ylim([-25 300] ) 
for n = 1:num_collects 
for m = l:num_lobes; 

% plot initial MUSIC DOAs 

line([repmat(coll_pos(1,n),1,q_est); endpt{n}(m,:,1)], ... 

[repmat{coll_pos{2,n),1,q_est);endpt{n}{m,:,2)],... 

'linestyle' , ' — ' ) 
if n > 1 && num_lobes ~=1; 

plot(intersect_pos{n}(m,:,1),intersect_pos{n}(m,:,2),... 

’o','color', colors(n,:)); 
if n > 2 

plot(mid_pos{n}(m,:,1),mid_pos{n}(m,:,2), 'ks' ); 

end 

end 

endpt_ff{n} (m,1:size{w{n},2), :) = w{n} (m, : , :)*300 + ... 

repmat{reshape(coll_pos(:,n),1,1,2) , [1,q_est,1]); 

%plot frequency corrected DOAs 

line([repmat{coll_pos{1,n),1,q_est); endpt_ff{n}(m,:,l)],... 
[repmat{coll_pos{2,n),1,q_est);endpt_ff{n}(m,:,2)]); 

end 

end 

if exist {' e_pos_est' ) 

plot{e_pos_est(1,:),e_pos_est(2,:), 'kd' , 'markersize' ,8) 

end 

xlabel{ 'In-track Range (nmi)' , 'fontsize' ,14) 
ylabel (' Cross-track Range (nmi) ', 'fontsize' ,14) 
title('Map View (frequecy fixed) ', 'fontsize ', 14) 
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% Transmitted signal plot 
fignum{4) = figure{); 
colors = get(gca, 'ColorOrder' ); 
for n = 1:q 

subplot(q,1,n) 

plot(t/le-9,real(sigs{n, :)), ’color',colors{n, ;)) 
if n < q 

set(gca, 'xtick’ , [ ]) 

end 

xlim([0 40]) 
legend(sig_names{n}) 

end 

xlabel{ 'time (ns)' , 'fontsize' ,14) 

ax = axes { ' Units ’ , ' Normal 'Visible’, ’off) ; 

set(get(ax, 'Ylabel '),' visibleon' ) 

ylabel( 'Volts (V)' , 'fontsize' ,14) 

title (' Transmitted signals’, 'fontsize’ ,14) 

set(get(ax, 'title ') , 'visible' , 'on' ) 

% Received signal plot 
fignum(6) = figureO; 
for n = 1:p 

subplot(p,1,n) 

plot(t/le-9,real(X(n, :))) 
if n < p 

set(gca, 'xtick’ , [ ]) 

end 

xlim([0 40]) 

legend(Strcat( 'r_’ ,num2str(n))) 

end 

xlabel( 'time (ns)' , 'fontsize' ,14) 

ax = axes ( ' Units ’ , ' Normal ’, 'Visible', ’off) ; 

set(get(ax, 'Ylabel' ), 'visible' , 'on' ) 

ylabel( 'Volts (V)' , 'fontsize' ,14) 

title( 'Received signals','fontsize', 14) 

set(get(ax, 'title ') , 'visible' , 'on' ) 

% FFT of received signal 
fignum(7) = figureO; 
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RSIG = fft(X(l,:)); 
stem(df/le6,abs(RSIG)) 

xlim([fmin fmax]/le6);set(gca, 'xtick’ ,{fmin:le8:fmax)/le6); 
grid on 

xlabel{ 'frequency (MHz)’font size', 14) 
ylabel( 'Voits (V)' , 'fontsize' ,14) 
title (' Received signal', 'fontsize' ,14) 

% Reconstructed signal plot 
fignum(8) = figure(); 
for n = 1:q_est 

subplot(q_est,1,n) 

plot(t/le-9,real(sig_est{1} (n, :)), 'color', colors (n, :)) 
if n < q 

set(gca, 'xtick' , [ ]) 

end 

xiim([0 40]) 

legend(Strcat{ 's_' ,num2str(n))) 

end 

xlabel( 'time (s)' , 'fontsize' ,14) 

ax = axes( 'Units' , 'Normal ', 'Visible', 'off') ; 

set(get(ax, 'Ylabel' ), 'visible' , 'on' ) 

ylabel( 'Volts (V)' , 'fontsize' ,14) 

set(get(ax, 'Ylabel' ), 'visible' , 'on' ) 

title (' Reconstructed signals', 'fontsize ', 14) 

set(get(ax, 'title ') , 'visible' , 'on' ) 

% Reconstructed Signal FFT plot 

fignum(9) = figureO; 

stem (df/le6, abs (SIG_est {1} . ' ) ) 

xlim([fmin fmax]/le6);set(gca, 'xtick' ,(fmin:le8:fmax)/le6); 
grid on 

xlabel( 'frequency (MHz) ', 'fontsize', 14) 

legend( 's_l' , 's_2' , ' s_3' , 's_4' ) 

ylabel( 'Volts (V)' , 'fontsize' ,14) 

title (' Reconstructed signals', 'fontsize ', 14) 

% Array beamshape plot at desired angle (0 for boresight) 

[~, indx] = min(abs(angles-O)) ; 
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CO = V{:,indx)’*v; 
figure {) 

plot(sinespace,20*logl0 (abs(CO))) 
ylim([-20 0]) 

title('Array Pattern', 'fontsize' ,14) 

xlabel{ 'sinespace' ) 

ylabel {' Normalized gain (dB) ') 


B.2 Signal Configuration Example Script 


o, g_ 
0 0 


o, 

0 


o 

0 


o, 
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o 
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o, 
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o, 
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o, 
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o, 
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o, 
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o, 

0 


o, 
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Signal Config Example 

This script defines receiver and signal parameters for use with the 
main MUSIC direction finding and geolocation script and provides 
example definitions. 

Receiver parameters (recParam.xxx): 

- centerfreq - array center frequency in Hz 

- sampfreq - sampling frequency in Hz (keep in mind top end of 
desired frequency range) 

- dur - signal duration (receiver duration will be slightly larger 
and will depend on array baseline) 

Signal parameters (sigParam{n}.xxx): 

- All signals require: 

- f - signal center frequency 

- type - options are 'sin', 'QPSK', 'LFM' 

- pos - location in the coordinate system. Should be positive y 

- snr - signal power (in dB); noise power is set to 0 dB if 
noise is enabled, so this effectively also sets the SNR 

- name - for plot labeling 

- LFM also requires: 

- B - signal bandwidth in Hz 

- QPSK also requires: 

- B - OPTIONAL, use if defining code as random sequence in order 
to obtain appropriate length 

- code - integer (0-3) based coding vector for QPSK phases. 
Bandwidth is based on code length: sigParam.B = K/recParam.dur, 


80 








29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 


where K is the number of subcodes in the vector. 

o, 

0 

% History: 

% DATE (YYMMDD) AUTHOR/EDITOR EMAIL NOTE 

% 141203 Chris Straessle gcstraes@nps.edu Original thesis code 

o, 
o 

OOOOOOQQQOOOOOOOOQQOOOOOOOOQQQOOOOOOOOQQQOOOOOOOOQQOOOOOOOOQQQOOOOO 


'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O 

% version name useful for specifying 
sigversion = 'sig_version' ; 

% receiver parameters 
recParam.centerfreq = le9; % Hz 
recParam.sampfreq = 4e9; % Hz 

recParam.dur = lOOe-9; % s 

% signal definitions 


sigParam) 

:i} 

.f = 

lOOOeG; 

sigParam) 

:i} 

.B = 

5 0 e 6 ; 

sigParam] 

:i} 

.type 

= 'QPSK'; 

sigParam] 

:i} 

. code 

= randi([0 3],floor 

sigParam] 

:i} 

. pos 

= [sind(15); cosd(15 

sigParam] 

:i} 

. snr 

= 10; 

sigParam] 

:i} 

. name 

= strcat(num2str(si 

sigParam) 

:2} 

.f = 

1105e6; 

sigParam) 

:2} 

.type 

= 'sin'; 

sigParam) 

:2} 

. pos 

= [90; 120]; 

sigParam) 

:2} 

. snr 

= 10; 

sigParam) 

:2} 

. name 

= strcat(num2str(si 

sigParam) 

:3} 

.f = 

1050e6; 

sigParam) 

:3} 

.B = 

100e6; 

sigParam) 

:3} 

.type 

= 'LEM' ; 

sigParam) 

:3} 

. pos 

= [sind(-32); cosd(- 

sigParam) 

:3} 

. snr 

= 10; 

sigParam) 

:3} 

. name 

= strcat(num2str((s 


, '-' , num2str{{sigParam{3}.f+si 


particular plots 


(sigParam{1}.B*recParam.dur),1); 

)]* 110 ; 

gParam{l}.f/le6), ' MHz QPSK'); 

gParam{2}.f/le6), ' MHz sine' ); 

32)]*150; 

igParam{3}.f-sigParam{3}.B/2)/le6)... 
Param{3}.B/2)/le6), ' MHz LEM'); 
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B.3 Interferometric Signal Generation Function 


function [sigs, X, e_pos, sig_names, doas, N, t] =... 

interferoinetry_sig_gen(sigParam, recParam, selected_sigs, ... 
coll_pos, coll_vel, rec_array, noiseon) 

o, 

o 

% Purpose; 

% Generate signals at receiver based on pre-defined signal parameters, 
% receiver array layout, collector position, and emitter position. 

% Useful for subspace-based direction finding algorithms. 

o 

0 

% Inputs: 

% - sigParam - structure of signal parameters (see example file) 

% - recParam - structure of receiver parameters (see example file) 

% - selected_sigs - vector where a specific subset of signals from 

% sigParam can be specified, e.g. [1 3] 

% - coll_pos - collector position matrix in [x;y], subsequent 

% positions in subsequent columns 

% - coll_vel - collector velocity in [x;y]. Assumed constant for all 

% collections. 

% - rec_array - linear receiver element location definition; should 

% specify receiver locations as distance along a line from the origin 

% - noiseon - equals 0 if additive noise is disabled and 1 if AWGN is 

% enabled 

o, 

0 

% Outputs: 

% - sigs - q X N vector of source signals (q is number of signals, N 

% is the number of samples 

% - X - p X N vector of received signals (p is number of receivers) 

% - e_pos - emitter positions in [x;y], subsequent emitter positions 

% in subsequent columns 

% - sig_names - cell array of signal names 

% - doas - signal DOAs in degrees, c x q matrix (c is the number 

% of collections) 

% - N - number of samples in collection 

% - t - time vector for collection (assumed to start at time zero) 

o 

O 

% Custom functions called: 
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% - gen_sin - generate phase shifted sinusoidal signal at receivers 

% - gen_qpsk - generate phase shifted QPSK signal at receivers 

% - gen_chirp - generate phase shifted LFM signal at receivers 

o, 

o 

% History: 

% DATE (YYMMDD) AUTHOR/EDITOR EMAIL NOTE 

% 141203 Chris Straessle gcstraes@nps.edu Original thesis code 

o, 

o 

9-9-9-9-2-S-9-9'9-9-9'9-9'9-S-2-S-9-9'9'9-9'9-9-9-S-2-9'9-9'9'9-9'9-9-2-S-9-9'9-9'9'9-9-9-9'9-S-2-9'9-9-9'9-9-9-9-2-S-9-9'9-9-9'9-9-9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


global c 

p = length(rec_array); % number of receivers 

q = length(selected_sigs); % number of signals 

%% Individula signal generation 


b = 0; 


for m = selected_sigs; % cycle through desired signals 

b = b + 1; % counter 

switch sigParam{m}.type 
case 'sin' 

[sig{b}, doas(b), ~] = gen_sin(sigParam{m}, ... 
recParam, coll_pos, coll_vel, rec_array); 
case 'QPSK' 


[sig{b}, doas(b), ~] 
recParam, coll_pos, 
case 'LFM' 


= gen_qpsk(sigParam{m}, . . . 
coll_vel, rec_array); 


[sig{b}, doas(b), ~] = gen_chirp(sigParam{m}, ... 
recParam, coll_pos, coll_vel, rec_array); 

end 

e_pos(:,b) = sigParam{m}.pos; % emitter positions 
sig_names{b} = sigParam{m}.name; % signal names 


% adjust signal power 

sig{b} = sig{b}*sqrt(10^{sigParam{m}.snr/10) ) ; 

end 


% number of samples; created to capture full signal length regardless of 
% DOA (includes additional time for 90 deg signal to traverse the array) 
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N = ceil { (recParam. dur+ (max {rec_array) 9rsind(90) ./c) ) *recParam. sampf req) ; 
t = (0:N-1)/recParam,sampfreq; % time vector 

%% Signal combination at array 

% initializations 
sigs = zeros(q,N); 

Xo = zeros(p,N); 

% create transmitted signal matrix and received signal matrix 
for m = 1:q 

sigs(m,1:size{sig{m},2)) = sig{m}(l,:); 

Xo{:,1:size(sig{m},2)) = Xo(:,1:size{sig{m},2)) + sig{m}; 

end 

Xa = Xo; % for future use (attenuation) 

%% Additive Noise 

% receiver noise variance (assume same for all elements) 
noise_var =1; % normalized noise 

if noiseon; 

rec_noise = sqrt (noise_var/2 ) 9c (randn (p, N) + 1 j 9rrandn (p, N) ) ; % AWGN 

else 

rec_noise = 0; % no noise case 

end 

X = Xa + rec_noise; % received signal matrix with receiver noise 


B.4 DOA and Time Delay Calculation Function 


function [doa, t_d] = doa_calc(sigParam, c_pos, c_vel, rec_array) 

o, 

o 

% Purpose: 

% Determine emitter DOA and delay vector which defines when the signal 
% (assumed plane wave) arrives at a receiver 

o, 

0 

% History: 

% DATE (YYMMDD) AUTHOR/EDITOR EMAIL NOTE 
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% 141203 Chris Straessle gcstraes@nps.edu Original thesis code 

o 

0 

9-9--9-9-9-9-9-9'9-9-9'2-9--9-9-9-9-9-9'9'9-9'9-9--9-9-9-9'9-9'9'9-9'9-9--9-9-9-9'9-9'9'9-9-9-9--9-9-9-9'9-9-9'9-9--9-9--9-9-9-9'9-9-9'Q-9--9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

global c 

% If no velocity, assume platform facing +x 
if norm{c_vel) == 0; 
c_vel = [ 1;0]; 

end 

e_pos = sigParam.pos; % emitter positions 
e_dist = norm(e_pos-c_pos); % distance to emitter 

% emitter DOA relative to array broadside 

doa = 90 - acosd({c_vel'*(e_pos-c_pos))./(e_dist*norm(c_vel))); 

% delay vector (column) 
t_d = rec_array*sind(doa)/c; 


B.5 Sinusoid Generation Function 


function [s, t_rec, doa, t_d] =... 

gen_sin(sigParam, recParam, c_pos, c_vel, rec_array) 

o, 

0 

% Purpose: 

% Generate sinusoidal pulse at receivers. 

o, 

o 

% Custom functions called: 

% - doa_calc - calculates signal DOA and delay vector 

o, 

o 

% History: 

% DATE (YYMMDD) AUTHOR/EDITOR EMAIL NOTE 

% 141203 Chris Straessle gcstraes@nps.edu Original thesis code 

o 

0 

9'9'9-S-9-S-9-9'9-9'9'9-9-9-S-2-9-9-9'9'9-9'9-9-9-9-9-9'9-9'9'9-9'9-9-2-S-2-9'9-9'9'9-9-9-9'2-S-2-9'9-9-9'9-9'9-9-9-S-9-9'9-9-9'9-9-9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

f = sigParam.f; % center frequency 
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ts = 1/recParam.sampfreq; % sample duration 
T = recParam.dur; % signal duration 

% calculate DOA and delay vector (column) 

[doa, t_d] = doa_calc(sigParam, c_pos, c_vel, rec_array); 

% collection time matrix, begins at time the first receiver receives 
% beginning of signal until the last receivers receives end of signal. 
% Signal at origin is time zero 

t_coll = repmat((ceil(min(t_d)/ts)*ts:ts:T-ts+max{t_d)),.. . 

[length(t_d),1]); 

% function to calculate delay at each each receiver 
t_ref = @(t)t - repmat(t_d,[1,length(t_coll)]); 

% matrix of shifted signal arrival times per receiver 
t_rec = t_ref(t_coll); 

% generate complex analytic sine signal 
s = exp (1 j* (2-^pi9rf + t_rec-pi/2) ) ; 

% set portions outside of pulse to zero 
s(t_rec < 0) = 0; 
s(t_rec > T-ts) = 0; 


B.6 QPSK Generation Function 


function [s, t_rec, doa, t_d] =... 

gen_qpsk(sigParam, recParam, c_pos, c_vel, rec_array) 

o, 

o 

% Purpose: 

% Generate QPSK pulse at receivers. 

o, 

o 

% Custom functions called: 

% - doa_calc - calculates signal DOA and delay vector 

o, 

0 

% History: 

% DATE (YYMMDD) AUTHOR/EDITOR EMAIL NOTE 

% 141203 Chris Straessle gcstraes@nps.edu Original thesis code 

o, 

0 
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9-9--9-S-9-S-9-9'9-9'9'9-9-9-S-2-S-9-9'9'9'9'9-9--9-S-9-9'9-9'9'9-9'9-9-9-S-2-9'9-9'9'9'9-9-9--2-S-2-9'9-9'9'9'9--9-9-2-S-9-9'9-9'9'9'9-9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


f = sigParam.f; % center frequency 

ts = 1/recParam.sampfreq; % sample duration 

T = recParam.dur; % signal duration 

z = sigParam.code; % code 

N = round{T/ts); % number of samples 

% determine appropriate subcode value for each sample 
z_long = repmat(z,[1 ceil(T/{length(z)*ts))]) 
z_long = z_long(:); 
z_long = z_long(l:N); 

% calculate DOA and delay vector (column) 

[doa, t_d] = doa_calc(sigParam, c_pos, c_vel, rec_array); 

% force signal to exist in positive time 
if min(t_d) < 0 

t_d = t_d - min(t_d); 

end 

% convert code integers to phase 
phi_opt = [pi/4, 3*pi/4, 5*pi/4, 7*pi/4]; 
phi = phi_opt(z_long+l); 

% collection time matrix, begins at time the first receiver receives 
% beginning of signal until the last receivers receives end of signal. 
% Signal at origin is time zero 

t_coll = repmat((ceil(min(t_d)/ts)*ts:ts:T-ts+max(t_d)),.. . 

[length(t_d),1]); 

% function to calculate delay at each each receiver 
t_ref = @(t)t - repmat(t_d,[1,length(t_coll)]); 

% matrix of shifted signal arrival times per receiver 
t_rec = t_ref(t_coll); 

% generate complex analytic signal 
s = exp (1 j* (2-^pi*f + t_rec) ) ; 

% set portions outside of pulse to zero 
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0; 


s(t_rec < 0) = 0; 
s(t_rec > T-ts) = 

% apply phase shifts to signal 
for m = 1:length(t_d) 

phaseshift = [zeros(1, abs(ceil(min{t_d)/ts)-ceil(t_d(m)/ts)))... 
exp(-1j*phi) . . . 

zeros(1,ceil(max(t_d)/ts)-ceil(t_d(m)/ts))]; 
s(m,:) = s(m,:).*phaseshift(1:length(s)); 

end 


B.7 LFM Generation Function 


function [s, t_rec, doa, t_d] =... 

gen_chirp(sigParam, recParam, c_pos, c_vel, rec_array) 

o 

O 

% Purpose; 

% Generate LFM pulse at receivers. 

o, 

o 

% Custom functions called: 

% - doa_calc - calculates signal DOA and delay vector 

o, 

0 

% History: 

% DATE (YYMMDD) AUTHOR/EDITOR EMAIL NOTE 

% 141203 Chris Straessle gcstraes@nps.edu Original thesis code 

o, 

0 

9'9'9-S-9-S-9-9'9-9'9'9-9-9-S-2-9-9-9'9'9'9'9-9'9-S-9-9'9-9'9'9'9'9-9-9-S-2-9'9-9'9'9-9-9-9'2-S-9-9'9-9'9'9'9'9-9-9-S-9-9'9-9-9'9'9-9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

f = sigParam.f; % center frequency 
ts = 1/recParam.sampfreq; % sample duration 
T = recParam.dur; % signal duration 
B = sigParam.B; % signal bandwidth 

% calculate DOA and delay vector (column) 

[doa, t_d] = doa_calc(sigParam, c_pos, c_vel, rec_array); 

% collection time matrix, begins at time the first receiver receives 
% beginning of signal until the last receivers receives end of signal. 
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% Signal at origin is time zero 

t_coll = repmat((ceil(min(t_d)/ts)*ts:ts:T-ts+max{t_d)),.. . 
[length(t_d),1]); 

% function to calculate delay at each each receiver 
t_ref = @(t)t - repmat(t_d,[1,length(t_coll)]); 

% matrix of shifted signal arrival times per receiver 
t_rec = t_ref(t_coll); 

% generate chirp phases 

phase = (f - B/2) .*(t_rec) + B/ (B^d) .9r(t_rec) .^2; 

% generate signal 
s = exp(1j*2*pi.*phase) ; 

% set portions outside of pulse to zero 
s(t_rec < 0) = 0; 
s(t_rec > T-ts) = 0; 


B.8 Subspace Calculation Function 


function [q_est, EigVec_s, EigVec_n, EigVal] = subspace_calc(Xf, ... 
sig_num_est, q, n) 

o, 

0 

% Purpose: 

% Calculate array covariance matrix (in frequency domain) and conduct 
% eigenvector decomposition of said matrix. Contains option to estimate 
% number of signals via Wax's MDL or Nadakuti and Edelman's NEW method, 

% or to force known number of signals 

o, 

o 

% Inputs: 

% - Xf - frequency domain version of received signal matrix 

% - sig_num_est - user parameter which defines how to estimate signal. 

% Options are 'none', 'MDL', or 'NEW' 

% “ q “ number of generated source signals 

% - n - collection number (used in alert message) 

o, 

0 

% Outputs: 

% - q_est - estimated number of source signals 
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% - EigVec_s - signal eigenvectors 

% - EigVec_n - noise eigenvectors 

% - EigVal - eigenvalues 

o, 

o 

% Custom functions called: 

% - num_sig_est_MDL - Wax's MDL estimator (1985) 

% - num_sig_est_NEW - Nadakuti and Edelman's estimator (2008) 

o, 

o 

% History: 

% DATE (YYMMDD) AUTHOR/EDITOR EMAIL NOTE 

% 141203 Chris Straessle gcstraes@nps.edu Original thesis code 

o 

0 

9-9'9-S-2-S-9-9'9-9-9'9-9-9-S-2-S-9-9'9'9-9'9-9-9-S-2-9'9-9'9'9'9'9-9-2-S-9-9'9-9'9'9-9-9-9'9-S-2-9'9-9-9'9-9-9-9-2-S-9-9'9-9-9'9-9-9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%% Subspace determination 
[p,N] = size(Xf) ; 

Rxx = Xf*Xf'/N; % Array covariance matrix 

[EigVec ,EigVal] = eig(Rxx); % eigendecomposition 

% Sort eigenvalues, then rearrange eigenvectors to match 
[EigVal,indx] = sort(diag(EigVal), 1, 'descend'); 

EigVec = EigVec ( :,indx); 

%% Signal Number Estimation 

switch sig_num_est 
case 'none' 

q_est = q; 
case 'MDL' 

% use Wax's MDL method 

q_est = num_sig_est_MDL(p, N, EigVal); 
case 'NEW' 

% use Nadakuti and Edelman's method 
q_est = num_sig_est_NEW(p, N, EigVal); 
otherwise 

error (' Incorrect signal number estimator specified!') 

end 
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if q_est q 

fprintf ([' Incorrect signal number estimate (%ld vice %ld)' ... 

' at collection number %ld.\n\n'], q_est, q, n); 
input {'Reset to known number of signals (y/n)?' ,reset_ans); 

switch reset_ans 

case { 'y' , 'Y' , 'yes' , 'Yes' , 'YES' } 

fprintf (' Resetting to known value\n\n') 
q_est = q; 
otherwise 

fprintf( 'Continuing without reset\n\n' ) 

end 

end 

%% Subspace determination cont. 

EigVec_s = EigVec(:,1:q_est); % signal eigenvectors 

EigVec_n = EigVec(:,q_est+l:p); % noise eigenvectors 


B.9 MUSIC Calculation Function 


function music_spec = inusic_calc (A, E_n) 

o 

0 

% Purpose: 

% Calculate the MUSIC (Schmidt, 1986) spectrum, given a steering 
% matrix A and the noise eigenvectors E_n of the received signal 
% covariance matrix. The function searches through provided steering 
% matrix to find steering vectors (i.e. angles) which are perpendicular 
% to the noise eigenvectors. These locations indicate signal directions 
% of arrival. 

o 

0 

% Inputs: 

% - A - steering matrix (n x m), where n is the number of receivers 

% and m is the number of discrete angles for which the spectrum is 

% being calculated. 

% - E_n - noise eigenvectors, (k x k), where 0 < k < n. 
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o, 

0 

% Outputs: 

% - music_spec - MUSIC spectrum 

o, 

o 

% History: 

% DATE (YYMMDD) AUTHOR/EDITOR EMAIL NOTE 

% I4I203 Chris Straessle gcstraes@nps.edu Original thesis code 

o, 

o 

9-9-9-9-2-S-9-9'9-9-9'9-9'9-S-2-S-9-9'9'9-9'9-9-9-S-2-9'9-9'9'9-9'9-9-2-S-9-9'9-9'9'9-9-9-9'9-S-2-9'9-9-9'9-9-9-9-2-S-9-9'9-9-9'9-9-9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

% execute MUSIC search 
for a = I:size(A,2) 

music_spec(I,a) = (A(:,a) '*A(:,a))/(A{:,a) '*(E_n*E_n')*A(: , a)) ; 

end 


B.IO 2D Interferometric Geolocation Function 


function e_pos_est = geolocate2(coll_pos, doa_est_ss) 

o, 

o 

% Purpose: 

% least squares solution to location based on closest point to multiple 
% DOA vectors. 

o, 

0 

% reference: 

% - http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html 
% NOTE: refer to thesis text for the 2D matrix version of formula in 
% MathWorld reference 

o, 
o 

% History: 

% DATE (YYMMDD) AUTHOR/EDITOR 
% 141203 Chris Straessle 

o, 
o 

9'9-9-S-2-S-9-9'9-9-9'9-9'9-S-2-S-9-9'9'9'9'9-9-9-S-2-9'9-9'9'9-9'9-9-2-S-2-9'9-9'9'9-9-9-9-2-S-2-9'9-9'9'9-9-9-9-2-S-9-9'9-9-9'9-9'9- 
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

[~,c] = size(coll_pos); % number of collections 

[~,q] = size(doa_est_ss); % number of emitters 

a = coll_pos; % collector position 


EMAIL NOTE 

gcstraes@nps.edu Original thesis code 
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% unit vector along DOA 
d = [permute(doa_est_ss,[3 1 2]); 

sqrt{1-(permute{doa_est_ss,[3 1 2])).^2)]; 

% create matrix problem for each signal (refer to thesis) 
for m = 1:q 

e = (d ( : , : , m) . ^2) ; 
f = d (1, : ,m) . *d (2, : ,in) ; 

A = [suin(e(2,:)), -l*sum (f); 

-l*suin(f), sum (e (1, : ) ) ] ; 

a_long = reshape(a, [1,2*c]) ; 

bl = reshape([e(2[l,2*c]); 

b2 = reshape(e(l,:)], [l,2*c]); 

B = [sum(a_long.*bl); sum{a_long.*b2)]; 

% solve matrix problem for estimated emitter location 
e_pos_est(:,m) = A\B; 

end 


B.ll Point Motion Metric Function 


function motion = pt_motion(positions) 

o 

0 

% Purpose: 

% Metric for estimating which DOA is true DOA and which is ambiguity. 

o, 

o 

% Inputs: 

% - positions - matrix of points (either DOA intersections or 
% DOA mid-positions. Format is: {coll #}( x|y|(z), sig #, lobe #) 

o, 

o 

% History: 

% DATE (YYMMDD) AUTHOR/EDITOR EMAIL NOTE 

% 141203 Chris Straessle gcstraes@nps.edu Original thesis code 

o 

0 

9-9'9-S-2-S-9-9'9-9'9'9-9-9-S-2-S-9-9'9'9'9'9-9'9-S-9-9'9-9'9'9'9'9-9-9-S-2-9'9-9'9'9'9-9-9'2-S-2-9'9-9'9'9'9'9-9-9-S-9-9'9-9'9'9'9-9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 
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% this allows metric to be used for either DOA intersections or 
% intersection mid-positions 
if isempty{positions{2}) 

first_pt =3; % mid-positions will start at collection 3 

else 

first_pt = 2 } % intersections will start at collection 2 

end 

% Change in x value from first to most recent positions 
x_diff = positions{ end} 1) - positions{first_pt1); 

% Change in y value from first to most recent positions 
y_diff = positions{ end} 2) - positions{first_pt}(:,:,2); 

% Magnitude (squared) of movement between first and last positions 
motion = x_diff.^2 + y_diff.'^2; 


B.12 Wax’s MDL Estimator Function 


function q_est = num_sig_est_MDL(p, N, EigVal) 

% Number of Signals Estimation (Wax 1985), MDL estimator 

o, 

o 

% p is the number of receivers 
% N is the number of samples 

% EigVal is the list of eigenvalues of the received signal covariance 
% matrix, sorted from largest to smallest. 

o, 

0 

% History: 

% DATE (YYMMDD) AUTHOR/EDITOR EMAIL NOTE 

% 141203 Chris Straessle gcstraes@nps.edu Original thesis code 

o, 

o 

9-9'9-9-9-9-9-9'9-9-9'9-9'9-9-9-9-9-9'9'Q-9'9-9'9-9-9-9'9-9'9'9-9'9-9'9-9-9-9'9-9'9'9-9'9-9'9-9-2-9'9-9-9'9-9'9-9'9-9-9-9'9-9-9'Q-9'9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

for k = 0:p-l; 

gO = prod(EigVal(k+1:p).^(1/(p-k))); 
aO = (1/(p-k)*sum(EigVal(k+1;p))) ; 
free_adj_param = .5*k*(2*p-k)*log(N); 

MDL(k+l) = -((p-k)*N)*log(gO/aO) + free_adj_param; 

end 
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[-f q_est] = min(MDL); 

q_est = q_est-l; % correct for MATLAB 1 vice 0 index 

end 


B.13 Nadakuti and Edelman’s Estimator Function 


function q_est = num_sig_est_NEW(p, N, EigVal) 

% Number of Signals Estimation (Nadakuditi and Edelman 2008), 'new' 

% estimator 

o, 

o 

% p is the number of receivers 
% N is the number of samples 

% EigVal is the list of eigenvalues of the received signal covariance 
% matrix, sorted from largest to smallest. 

o 

0 

% History: 

% DATE (YYMMDD) AUTHOR/EDITOR EMAIL NOTE 

% 141203 Chris Straessle gcstraes@nps.edu Original thesis code 

o, 

o 

9 - 9 ' 9 - 9 - 9 - 9 - 9 - 9 ' 9 - 9 - 9 ' 2 - 9 ' 9 - 9 - 9 - 9 - 9 - 9 ' 9 ' 9 - 9 ' 9 - 9 ' 9 - 9 - 9 - 9 ' 9 - 9 ' 9 ' 9 - 9 ' 9 - 9 ' 9 - 9 - 9 - 9 ' 9 - 9 ' 9 ' 9 - 9 ' 9 - 9 ' 9 - 9 - 9 - 9 ' 9 - 9 - 9 ' 9 - 9 ' 9 - 9 ' 9 - 9 - 9 - 9 ' 9 - 9 ' 9 ' 9 - 9 ' 9 - 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

for k = 0:p-l; 

parti = (p-k) *sum (EigVal (k+l: p) . ^2)/sum (EigVal (k+l: p) ).''2 ; 
t_k = (parti - (l+p/N))*p; 

NEW(k+l) = .5*(N*t_k/p)^2 + 2*(k+l); 

end 

[~, q_est] = min(NEW); 

q_est = q_est-l; % correct for MATLAB 1 vice 0 index 

end 
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